REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments 

regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 

Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA,  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  oenalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1.  REPORT  DATE  (DD-MM-YYYY) 
29-09-2012 


4.  TITLE  AND  SUBTITLE 

Detection  of  Low-order  Curves  in  Images  using 
Biologically-plausible  Hardware 


2.  REPORT  TYPE 
Final  Report 


3.  DATES  COVERED  (From  -  To) 

1 -Oct-201 1  -  30-Jun-2012 


5a.  CONTRACT  NUMBER 
W911NF-11-1-0271 


5b.  GRANT  NUMBER 


6.  AUTHORS 
Wesley  E.  Snyder 


5c.  PROGRAM  ELEMENT  NUMBER 
611102 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

North  Carolina  State  University 
Research  Administration 
NC  State  University 

Raleigh,  NC  27695  -7514 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND 
ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

60356-NS-II.l 


12.  DISTRIBUTION  AVAILIBILITY  STATEMENT 
Approved  for  Public  Release;  Distribution  Unlimited 


13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


14.  ABSTRACT 

Research  was  conducted  to  explore  the  possibility  that  within  the  human  visual  cortex  there  lies  specialized 
networks  which  function  in  a  manner  similar  to  that  of  the  Hough  transform.  Simulation  software  was  written  and 
tested  on  images.  The  tests  were  compared  by  results  reported  by  psychologists  testing  human  performance,  and 
found  to  produce  consistent  results.  The  conclusion  is  made  that  there  is  no  reason  to  rule  out  the  possibility  of 
Hough-like  circuits  in  the  brain  which  detect  straight  lines  in  the  peripheral  visual  field. 


15.  SUBJECT  TERMS 

Neural  modeling,  Hough,  straight  line  detection 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

15.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

ABSTRACT 

OF  PAGES 

Wesley  Snyder 

UU 

UU 

UU 

UU 

19b.  TELEPHONE  NUMBER 

919-515-5114 

Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39. 1 8 


Report  Title 

Detection  of  Low-order  Curves  in  Images  using  Biologically-plausible  Hardware 

ABSTRACT 

Research  was  conducted  to  explore  the  possibility  that  within  the  human  visual  cortex  there  lies  specialized  networks  which  function  in  a 
manner  similar  to  that  of  the  Hough  transform.  Simulation  software  was  written  and  tested  on  images.  The  tests  were  compared  by  results 
reported  by  psychologists  testing  human  performance,  and  found  to  produce  consistent  results.  The  conclusion  is  made  that  there  is  no 
reason  to  rule  out  the  possibility  of  Hough-like  circuits  in  the  brain  which  detect  straight  lines  in  the  peripheral  visual  field. 


Enter  List  of  papers  submitted  or  published  that  acknowledge  ARO  support  from  the  start  of 
the  project  to  the  date  of  this  printing.  List  the  papers,  including  journal  references,  in  the 
following  categories: 

(a)  Papers  published  in  peer-reviewed  journals  (N/A  for  none) 


Received  Paper 

TOTAL: 

Number  of  Papers  published  in  peer-reviewed  journals: 


(b)  Papers  published  in  non-peer-reviewed  journals  (N/A  for  none) 


Received  Paper 


TOTAL: 

Number  of  Papers  published  in  non  peer-reviewed  journals: 


(c)  Presentations 


Number  of  Presentations:  0.00 

Non  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 

Received  Paper 

TOTAL: 

Number  of  Non  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 

Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 

Received  Paper 

TOTAL: 

Number  of  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


(d)  Manuscripts 


Received 


TOTAL: 

Number  of  Manuscripts: 

Books 

Received  Paper 

TOTAL: 

Patents  Submitted 


Patents  Awarded 

Awards 


Graduate  Students 

PERCENT  SUPPORTED  Discipline 
0.47 

0.47 

1 

Names  of  Post  Doctorates 


Names  of  Faculty  Supported 


NAME 

PERCENT  SUPPORTED  National  Academy  Member 

Wesley  Snyder 

0.20 

FTE  Equivalent: 

0.20 

Total  Number: 

1 

Names  of  Under  Graduate  students  supported 


NAME 

PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

NAME 
Theju  Jacob 

FTE  Equivalent: 
Total  Number: 


Student  Metrics 

This  section  only  applies  to  graduating  undergraduates  supported  by  this  agreement  in  this  reporting  period 

The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period: .  0.00 

The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period  with  a  degree  in 

science,  mathematics,  engineering,  or  technology  fields: .  0-00 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  continue 

to  pursue  a  graduate  or  Ph.D.  degree  in  science,  mathematics,  engineering,  or  technology  fields: 0  00 

Number  of  graduating  undergraduates  who  achieved  a  3.5  GPA  to  4.0  (4.0  max  scale): .  0.00 

Number  of  graduating  undergraduates  funded  by  a  DoD  funded  Center  of  Excellence  grant  for 

Education,  Research  and  Engineering: .  o.OO 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  intend  to 

work  for  the  Department  of  Defense .  0.00 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  receive 

scholarships  or  fellowships  for  further  studies  in  science,  mathematics,  engineering  or  technology  fields: .  0.00 

Names  of  Personnel  receiving  masters  degrees 

NAME 

Total  Number: 

Names  of  personnel  receiving  PHDs 

NAME 

Total  Number: 

Names  of  other  research  staff 

NAME  PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

Sub  Contractors  (DD882) 


Inventions  (DD882) 


See  Attachment. 


Scientific  Progress 
Technology  Transfer 


Detection  of  Low-order  Curves  in  Images  using  Biologically-plausible  Hardware 


FINAL  REPORT 

Submitted  to  the  Army  Research  Office 
by 

North  Carolina  State  University 


Wesley  Snyder 
Principal  Investigator 


Covering  the  period  1  October,  2011  -  1  July  2012 


1 


Contents 


1  Project  Description  5 

2  Objectives  6 

3  Technical  Background  8 

3.1  The  Standard  Model .  8 

3.2  The  log-polar  map .  9 

4  An  Architecture  for  Non-local  Vision  11 

4.1  Background:  the  Hough  Transform .  13 

4.2  The  Architecture .  16 

5  Approach  18 

5.1  Estimating  Orientations .  18 

5.1.1  Estimating  Orientation  Using  the  Traditional  Hough  Transform  (mode  0)  .  .  19 

5.1.2  Estimating  Orientation  By  Interpolation  (mode  1) .  19 

5.1.3  Estimating  Orientation  Using  the  Pseudo-inverse  (mode  2)  19 

5.1.4  Estimating  Orientation  Using  the  maximum  Gabor  Output  (mode  3) .  20 

5.1.5  Estimating  Orientation  by  Incrementing  all  the  Directions  (mode  4) .  20 

5.1.6  Estimating  Orientation  by  Fitting  a  Parabola  (mode  5)  .  20 

5.1.7  Estimating  Orientation  by  Linearly  Interpolating  the  Strongest  Two  Direc¬ 
tions  (mode  6)  .  20 

5.2  Accumulation .  21 

6  Results  26 

6.1  The  Experiments .  26 

6.2  Detection  as  a  Function  of  Length  .  27 

6.3  Detection  as  a  Function  of  Microline  Orientation .  27 

6.4  Detection  as  a  Function  of  Clutter .  34 

6.5  Comparison  with  human  experiments .  34 

6.5.1  Similarities  in  experimental  parameters  .  37 

6.5.2  Differences  in  experimental  parameters .  37 

6.5.3  Results  Compared  with  Human  Performance  in  Literature .  37 


2 


7  Future  Work  41 

7.1  Correlation  with  Human  Behavior  .  41 

7.2  Hardware .  41 

7.3  Learning .  42 

7.4  Impact  on  Models  for  the  Visual  Cortex .  42 

8  An  Algorithm  for  Drawing  Straight  Lines  44 

9  Conclusion  46 

A  Determining  Orientation  from  Four  Estimates  47 

B  The  Log-Polar  Transform  48 

B.l  Introduction .  48 

B.2  Beginnings  -  1960s,  1970s .  48 

B.3  1980s, 1990s, 2000s .  51 

B.3.1  Emphasis  of  Central  Vision  in  Retina  .  53 

B.3. 2  Emphasis  of  Central  Vision  in  Retino  Cortical  Pathway .  54 

B. 3. 3  Later  studies,  use  of  fMRI .  56 

B. 4  Conclusion  .  60 

C  Software  61 

C. l  Definitions  and  structures .  63 

C.2  Help  Function .  65 

C.3  Parsing  the  Command  Line .  66 

C.4  Saving  Accumulator  to  Disk .  68 

C.5  Computing  a  histogram  of  the  Accumulator  values .  69 

C.6  Main  Program  .  70 

C.7  saveacc .  74 

C.8  Find  peaks  in  accumulator .  76 

C.9  printmatrix .  79 

C.10  Initialization  .  80 

C.ll  Set  up  connections  of  Axons .  86 

C.12  High  Pass  Filtering .  87 

C. 12.1  Dynamic  Response  of  the  input  system  .  87 

C.13  Non  Maximum  Suppression .  89 

C.14  Calling  the  Gabor  Edge  Detectors  .  91 

C.l 5  Simulation  of  the  Accumulator .  94 

C.16  incacc . 102 

C.l 7  Low-pass  filtering  the  Accumulator . 104 

C.18  Making  a  Gaussian  blurring  kernel . 107 

C.l 9  Set  an  image  to  all  zeros . 108 

C.20  Initializing  an  Accumulator . 109 


3 


C.21  Disk  reread . 110 

C.22  Writing  a  file . Ill 

C.23  Test  for  image  null . 112 

C.24  Clip  a  floating  point  image  . 113 

C.25  Blurring  the  vector  of  angle  measurements . 114 

C.26  Display  the  Inverse  Accumulator . 115 

C.27  MakeFake  Function  only  for  Testing . 116 

C.28  Accumulator  Operations  Functions . 116 

C.28.1  Video  Scale  Accumulator . 116 

C.28. 2  Mark  POints  of  Interest  in  Accumulator . 117 

C.28. 3  ShowAccNeighborhood . 117 

C.28. 4  Determine  if  a  point  is  a  Local  Maximum . 118 

C.28. 5  Make  an  Accumulator . 118 

C.28. 6  test  Indexing . 119 

C.28. 7  Cosine  of  an  Angle  Specified  in  Degrees . 120 

C.28. 8  Since  of  an  Angle  Specified  in  Degrees . 120 

C.28. 9  Convert  an  Accumulator  to  an  IFS  Image . 120 

C.29  Reconstruct  and  Display  the  Original  Image . 124 


4 


Chapter  1 

Project  Description 


The  work  is  motivated  by  the  following  experience,  related  to  the  PI: 

I  was  looking  over  a  large  body  of  water  with  land  on  the  other  side.  The  sky  was 
cloudless  and  blue.  Suddenly,  I  was  aware  of  a  flash  of  a  straight  line  in  the  upper 
right  of  my  visual  field.  When  I  shifted  my  gaze  to  that  point,  I  realized  a  bird  had 
momentarily  flown  in  such  a  way  as  to  create  a  straight  line  between  a  water  tower 
and  a  tree. 

This  simple  experience  illustrates  several  aspects  of  the  perception  of  straight  lines.  First,  collinear- 
ity  can  be  detected  in  non-foveal  areas  of  the  visual  field;  the  observer  was  not  looking  at  the  bird 
until  after  his  attention  was  shifted.  Second,  since  directed  attention  was  not  on  the  area  of  the 
visual  field  containing  the  line,  this  suggests  an  operation  performed  by  “hardware”  in  the  visual 
cortex.  Third,  this  is  an  unusual  phenomenon,  it  was  attention-diverting  because  the  background 
was  so  extremely  uncluttered,  arguing  for  a  system  that  is  sensitive  to  background  clutter. 

In  later  sections  of  this  report,  we  will  refer  back  to  this  experience. 

In  this  work,  we  have  sought  a  biologically-plausible  computing  architecture  which  can  identify 
simple  features  such  as  straight  lines  in  images.  We  sought  an  architecture  which  will  detect  these 
features  over  very  wide  expanses  of  the  visual  field,  not  restricted  to  the  fovea  or  a  local  receptive 
field.  We  sought  answers  to  the  following  questions: 

•  whether  mammals  have  “hardware”  for  straight  line  detection, 

•  what  the  nature  of  such  hardware  is,  given  that  only  low  precision  calculations  are  available, 

•  how  to  detect  and  quantify  such  functions  in  the  visual  cortex,  and 

•  how  to  build  an  electronic  real-time  version. 

We  have  concluded  that  a  computing  architecture  based  on  straightforward  image  analysis  opera¬ 
tions  provides  Straight  Line  Detection  (SLD)  consistent  with  human  behavior. 
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Chapter  2 

Objectives 


In  chapter  4.2  we  will  present  a  neural  architecture  which  is  reasonable  to  describe  the  high-speed 
SLD  experience  described  above.  We  have  simulated  this  architecture,  and  tested  its  performance. 
We  will  then  describe  future  work  performing  comparisons  with  carefully  conducted  human  behavior 
experiments  to  confirm  that  our  proposed  and  simulated  architecture  is  a  reasonable  explanation 
for  how  human  detect  straight  lines. 
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Figure  2.1:  A  Gabor  filter  combines  derivative  of  a  Gaussian  with  sinusoidal  modulation.  The 
sensitivity  of  the  edge  detection  neurons  is  consistent  with  a  Gabor  filter.  The  filter  shown  in  this 
figure  produces  an  estimate  of  the  first  derivative  of  a  brightness  edge. 
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Chapter  3 


Technical  Background 


Since  the  work  of  Hubei  and  Wiesel,  (now  collected  into  a  book  [38]),  it  has  been  possible  to 
represent  the  lower  levels  of  the  biological  vision  system  by  a  convolution  with  a  collection  of  edge- 
sensitive  kernels.  For  example,  we  know  by  experiment  that  corresponding  to  every  point  in  the 
retina,  there  are  cells  which  are  sensitive  to  edges/lines  with  varying  orientation. 

The  output  of  these  cells  may  be  simulated  by  convolution  with  an  oriented  kernel  such  as  a  Gabor 
filter.  These  feature  detectors  are  known  as  Simple  cells  at  layer  1  or  just  SI  cells. 


3.1  The  Standard  Model 


The  standard  model  for  vision  has  been  developed  in  a  variety  of  versions  over  a  number  of  years 
by  Tomaso  Poggio  and  colleagues.  This  particular  model  for  vision  has  always  been  very  careful 
to  avoid  using  any  architectural  components  that  could  not  be  found  in  a  biological  brain,  and  to 
substantiate  any  models  with  actual  neural  recordings  in  so  far  as  possible.  We  follow  the  same 
philosophy  in  this  work. 

A  thorough  explanation  of  the  Standard  Model  is  provided  in  [57].  The  model  begins  with  area  VI 
of  the  visual  cortex.  From  many  neurophysiological  experiments,  VI  is  known  to  be  mapped  in  a 
retinotropic  way  (neurons  close  to  each  other  in  VI  correspond  to  areas  in  the  visual  field  which 
are  likewise  close  to  each  other).  The  details  of  processing  done  by  ganglion  cells  in  the  retina,  the 
optic  nerve,  and  the  lateral  geniculate  nucleus  (LGN)  are  lumped  into  the  observation  that  there 
are  cells  in  VI  which  are  sensitive  to  edges  in  the  image.  Different  cells  are  sensitive  to  edges  with 
different  orientations.  There  seem  to  be  about  4  preferred  orientations,  horizontal,  vertical,  and  the 
two  45-degree  slanted  edges.  The  output  of  these  cells  is  consistent  with  convolution  of  the  image 
with  a  Gabor  filter  (See  Figure  2.1)with  that  specific  orientation.  The  convolution  operation  may 
be  implemented  in  simulation  by 


9  =  9 


(  T.Uv’rf  \ 


(3.1) 


where  the  wj  are  the  weights  corresponding  to  the  filter,  the  Xj  are  the  inputs  from  the  LGN,  and 
g  is  the  sigmoid  function  commonly  associated  with  neural  response  functions.  Refer  to  [57]  for 
details. 

The  outputs  of  the  SI  cells  are  inputs  to  Cl  cells,  which  have  larger  receptive  fields,  and  which  per¬ 
form  a  maximum  operation  over  a  local  neighborhood.  By  selecting  the  output  from  the  strongest 
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Figure  3.1:  LEFT  When  a  straight  or  gently  curving  line  is  occluded,  the  human  visual  system 
infers  the  existence  of  the  occluded  portions  of  the  line.  This  occurs  even  when  the  space  between 
the  collinear  segments  is  significantly  larger  than  the  observation  area  of  the  fovea.  RIGHT: A  line 
in  the  image  is  parameterized  by  its  distance  from  the  origin  (p),  and  the  angle  it  makes  with  the 
x  axis  (9).  Any  point  on  that  line  may  be  described  using  polar  coordinates  r  and  a. 


SI  cell,  and  effectively  suppressing  the  others,  a  small  degree  of  translational  invariance  is  provided, 
coupled  with  a  more  accurate  positioning  of  actual  edges.  This  is  similar  to  the  function  of  the 
Canny[17]  edge  operator,  which  uses  nonmaximum  suppression  to  locate  the  maximum  of  local 
edge  detector  responses.  Biologically-plausible  implementations  of  the  maximum  operation  may  be 
described  by  the  “softmax”  operation  [57]  satisfying  the  equation 


V  =  9 


\*+T,u*V  ' 


(3.2) 


The  maximum  operation  is  also  computable  using  a  locally-connected  network  of  neurons  using 
inhibition  to  implement  winner-take-all. 

Thus,  after  the  first  SI /Cl  pair,  edges  with  specific  orientations  are  identified  with  considerable 
precision.  This  will  become  the  input  to  our  SLD  network.  In  [57] ,  additional  levels  of  S-  and  C-like 
cells  are  proposed  at  higher  levels  of  the  brain,  with  increasing  size  and  generality,  and  shown  to 
provide  a  feedforward  network  that  can  explain  foveal  shape  recognition.  However,  the  SI- Cl  level 
is  all  that  is  required  for  the  SLD  operation  described  in  this  work. 


3.2  The  log-polar  map 

In  1977,  Eric  Schwartz  [55]  published  his  seminal  paper  illustrating  the  “log-polar  transform” .  This 
transform,  denoted  L  :  C  — »  C  describes  the  image  plane  by  an  isomorphism  to  the  complex  plane, 
so  that  a  point  with  coordinates  (x,  y)  in  the  image  is  represented  by  the  single  complex  variable 
2.  Similarly,  a  point  in  the  cortical  plane  is  represented  by  w  E  C,  and  the  log  polar  mapping  is 

w  =  ln(z)  (3.3) 
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This  representation  is  shown  to  be  consistent  with  a  substantial  body  of  neurophysiological  liter¬ 
ature,  and  several  advantages  result  from  this  representation.  For  example,  scale  change  in  the 
image  is  converted  to  a  simple  translation  on  the  cortex. 

The  complex  number  w  in  magnitude,  direction  form  is  denoted  w  =  (p,  6),  where 

p  =  logi/x2  +  y 2 
8  =  arctan  - 

X 


Figure  3.2:  Log  Polar  Transformation 


p  is  hence  the  logarithm  of  the  distance  to  a  given  point  from  the  origin  and  9  is  the  angle  between 
the  axis  and  the  line  joining  the  given  point  and  the  origin.  In  image  processing  applications,  the 
mapping  aids  in  data  reduction  by  reducing  the  resolution  at  image  boundaries.  The  transform  is 
also  widely  used  in  applications  like  image  registration,  tracking  and  target  recognition. 

In  Appendix  B,  the  transform  is  explained  in  more  detail.  This  analysis  was  provided  by  the 
graduate  student  supported  on  this  grant,  but  independent  of  the  grant  itself. 


10 


Chapter  4 


An  Architecture  for  Non-local  Vision 


A  line,  or  collection  of  nearly  collinear  line  segments,  in  an  image  may  be  detected  by  the  well-known 
Hough  transform  (HT)  [36]. 

The  HT  is  a  type  of  accumulator  method ,  and  since  the  original  paper,  many  variations  and  ex¬ 
tensions  have  been  made  [7,  19],  which  follow  the  Hough  philosophy  of  using  accumulation  to 
determine  consistency.  Some  of  these  methods  are  based  on  neural  models  [12,  13,  8,  9,  10,  11].  In 
this  presentation,  we  discuss  only  the  original  straight-line  version. 

Accumulator  methods  allow  consistent  but  spatially  diverse  image  segments  to  “vote”  for  consistent 
interpretations.  Because  the  voting  is  an  addition,  zero  mean  noise  is  averaged  out.  The  HT 
methods  have  been  shown  [33,  16]  to  provide  highly  robust  ways  to  combine  local  measurements 
to  make  global  decisions. 

1.  The  method  collects  input  from  spatially  diverse  features  in  the  image  and  increments  a 
single  point/neighborhood  in  the  transform  image1.  Thus,  non-connected  but  almost-collinear 
features  in  the  image  all  enhance  the  same  point  in  the  transform  space.  This  produces  a 
means  for  dealing  simply  with  occlusion.  Figure  3.1  shows  an  example  in  which  the  human 
extrapolates,  and  infers  the  presence  of  straight  lines. 

2.  The  accumulation  process  is  simply  additions.  The  method  is  therefore  robust  to  most  addi¬ 
tive  noise,  since  positive  and  negative  noise  cancels  out. 

3.  Features  which  are  identical  (e.g.  straight  lines)  but  located  in  different  points  in  the  image 
transform  to  different  points  in  the  transform  space. 

We  have  observed  that  the  advantages  of  accumulator-based  methods  are  quite  significant,  and  in 
this  work  we  insisted  on  finding  methods  to  use  them  in  a  biologically-plausible  network,  though 
there  may  be  alternatives.  For  example,  Basak  [8]  develops  a  method  for  finding  straight  lines  in 
images,  however,  no  accumulators  are  used  (although  the  title  of  the  paper  includes  the  words  Hough 
Transform).  Instead,  a  search  is  made  for  a  collection  of  weight  vectors,  each  vector  representing 
a  straight  line.  Neural  nets  methods  (weight  vector  estimation)  are  used  in  an  iterative  algorithm 
which  requires  that  the  maximum  allowable  number  of  vectors  be  known  a-priori.  In  contrast,  in 
this  work,  we  made  extensive  use  of  accumulators. 

Here,  we  consider  only  straight  lines,  although  we  have  done  considerable  research  in  the  use 
of  accumulator-based  techniques  in  detecting  general  shapes [41,  42],  In  Figure  3.1(RIGHT)  is 

1  Since  the  transform  is  implemented  using  a  2-D  array,  we  will  use  the  word  image. 
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Figure  4.1:  Each  edge-sensitive  cell  along  the  line  (p,0)  extends  a  single  connection  to  the  the 
accumulator  cell  characterizing  the  edge.  In  this  figure,  two  edges  are  illustrated,  with  orientations 
of  9\  and  62 .  The  accumulator  itself  is  organized  in  polar  coordinates,  so  the  angle  from  the  center 
of  the  accumulator  (O)  is  parallel  to  the  line  which  it  represents. 


12 


Figure  4.2:  An  image  which  is  the  output  from  an  edge  detector.  The  human  can  immediately 
discern  that  this  boundary  consists  of  two  straight  segments. 


illustrated  a  straight  line  parameterized  by  its  distance  from  the  origin  and  the  angle,  </>,  it  makes 
with  the  horizontal  axis. 


4.1  Background:  the  Hough  Transform 

Suppose  one  is  tasked  with  the  problem  of  finding  the  straight  lines  in  the  image  shown  in  Figure 
4.2.  If  only  one  straight  line  were  present  in  the  image,  we  could  use  straight  line  fitting  to  determine 
the  parameters  of  the  curve. 

We  could  represent  lines  by  equations  of  the  form 

y  =  ax  +  b  (4.1) 

But  we  have  two  line  segments  here.  If  we  could  segment  this  first,  then  we  could  fit  each  segment 
separately  -  yes,  this  is  a  segmentation  problem  -  but  we  are  segmenting  a  boundary  into  boundary 
segments  rather  than  segmenting  an  image  into  regions.  In  this  section,  we  will  learn  how  to  do 
this. 

First,  let  us  prove  an  illustrative  theorem. 

Definition :  given  a  point  in  a  d-space,  and  a  parametrized  expression  defining  a  curve  in  that 
space,  the  parametric  transform  of  that  point  is  the  curve  which  results  from  treating  the  point 
as  a  constant  and  the  parameters  as  variables.  For  example,  Eq.  4.1  produces  the  parametric 
transform 

b  =  y  —  xa  (4.2) 
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which  is  itself  a  straight  line  in  the  2-space  <  a,  b  >.  Given  the  point  x  =  3,y  =  5,  then  the 
parametric  transform  is  b  =  5  —  3a. 

Theorem:  If  n  points  in  a  2-space  are  collinear,  all  the  parametric  transforms  corresponding  to 
those  points,  using  the  form  b  =  y  —  xa  intersect  at  a  common  point  in  the  space  <  a,b  >. 

Proof:  Suppose  n  points  {xi,  y\),  (x2, 2/2), , ,  (xn,  yn}  all  satisfy  the  same  equation 

y  =  a0x  +  b0  (4.3) 

Consider  two  of  those  points,  ( Xi,yi )  and  ( Xj,yj )  .  The  parametric  transforms  of  the  points 
are  the  curves  (which  happen  to  be  straight  lines) 

yt  =  Xid  +  b  (4.4) 

y.j  =  Xja  +  b  (4.5) 

(4.6) 


which  we  rewrite  to  make  clear  the  fact  that  a  and  b  are  independent  variables 


b  =  yi  -  Xia  (4.7) 

b  =  yj  —  Xja  (4.8) 


The  intersection  of  those  two  curves,  straight  lines  in  a,  b  is  a  single  point. 

Solving  the  two  equations  of  Eq.  4.8  simultaneously  results  in 

yi  -  yj  =  (xi  -  Xj)a  (4.9) 


and  therefore  a  =  ■ 

We  substitute  a  into  Eq.  4.6  to  find  b, 


b  =  yt-  Xi 


Vi  ~  Vj 

Xi  —  Xj 


(4.10) 


and  we  have  the  a  and  b  values  where  the  two  curves  intersect.  However,  we  also  know  from  Eq. 
4.3  that  all  the  xs  and  ys  satisfy  the  same  curve.  By  performing  that  substitution  into  Eq.  4.10, 
we  obtain 


,  ,  ,  7  x  ((aoxj  -  60)  -  ( aoXi  +  b0)) 

b  =  ( a0Xi  +  b0)  -  Xi - - - — — — - 

3C  2  Xj 


(4.11) 


which  simplifies  to 


b  =  ( a0Xi  +  b0)  -  Xia0  =  b0 


(4.12) 


Similarly, 


Vi  ~  Vj  _  («o Xi  +  60)  -  ( a0Xj  +  b0) 


(4.13) 


Thus,  for  any  two  points  along  the  straight  line  parameterized  by  ao  and  bo,  their  parametric 
transforms  intersect  at  the  point  a  =  ao  and  b  =  bo-  Since  the  transforms  of  any  two  points 
transforms  intersect  at  that  one  point,  all  such  transforms  intersect  at  that  common  point.  QED. 
Review  of  concept:  Each  POINT  in  the  image  produces  a  CURVE  (possibly  straight)  in  the  pa¬ 
rameter  space.  If  the  points  all  lie  on  a  straight  line  in  the  image,  the  corresponding  curves  will 
intersect  at  a  common  point  in  parameter  space. 
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The  Problem  with  Vertical  Lines  In  order  to  deal  with  vertical  lines,  in  which  the  slope,  a,  is  infinite, 
we  choose  a  different  form  for  the  straight  line 

p  =  x  cos  0  -\-  y  sin  6  (4.14) 

For  a  given  pair  of  values  for  p  and  8,  the  set  of  points  which  satisfy  Eq.  4.14  can  be  shown  to  be 
a  straight  line. 

This  representation  of  a  straight  line  has  a  number  of  advantages.  Unlike  the  use  of  the  slope,  both 
of  these  parameters  are  bounded;  p  can  be  no  larger  than  the  largest  diagonal  of  the  image,  and 
8  need  be  no  larger  than  2n.  A  line  at  any  angle  may  be  represented  without  singularity.  The 
use  of  this  parameterization  of  a  straight  line  solves  one  of  the  problems  which  confronts  us,  the 
possibility  of  infinite  slopes.  The  other  problem  is  the  calculation  of  intersections. 

How  to  Find  Intersections  -  Accumulator  Arrays 

It  is  not  feasible  to  find  all  intersections  of  all  curves,  and  to  then  determine  which  of  those  are  close 
together.  Instead,  we  make  use  of  the  concept  of  an  accumulator  array.  To  create  an  accumulator 
array,  we  make  an  image,  say  360  columns  by  512  rows.  We  initialize  each  of  the  pixels  to  zero. 
From  now  on,  we  will  refer  to  the  pixels  of  this  special  image  as  accumulators.  Figure  4.4  illustrates 
plotting  two  straight  lines  through  an  accumulator  array  using  the  following  algorithm: 

•  For  each  point  (, Xi,yi )  in  the  edge  image,  do 

1.  for  all  values  of  8  compute  p. 

2.  at  the  point  p,  8  in  the  accumulator  array,  increase  the  value  at  that  point  by  an  amount 
proportional  to  the  strength  of  the  edge. 

This  algorithm  results  in  multiple  increments  of  those  accumulators  corresponding  to  intersections 
being  increased  more  often.  Thus  the  peaks  in  the  accumulator  array  correspond  to  multiple 
intersections,  and  hence  to  proper  parameter  choices. 


Figure  4.3:  Two  very  noisy  lines. 
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In  this  work,  we  do  not  use  the  Hough  trans¬ 
form  as  described  in  this  subsection  (except  for 
test  validation),  however,  we  do  make  use  of  ac¬ 
cumulators  to  find  consistent  contributions  to 
collinear  segments. 

4.2  The  Architecture 

There  are  two  arrays  of  interest,  the  image 
and  the  accumulator.  The  image  is  defined  by 
the  output  of  Cl  cells,  and  denoted  Cl(x,y,  9). 

The  spatial  coordinates  of  the  Cl  cell  could  be 
Cartesian  (x,  y )  or  polar  (r,  </>).  In  the  proposed  system,  either  may  be  used,  but  we  choose  Cartesian 
simply  because  it  is  easier  for  humans  to  visualize.  We  choose  to  use  the  output  of  Cl  cells  because 
those  cells  reflect  the  results  of  two  lower-level  stages  of  processing,  convolution  and  softmax.  The 
convolution  is  performed  by  SI  cells  as  described  in  section  3.1,  and  similarly,  the  Cl  cells  provide 
a  degree  of  invariance  to  small  rotations  and  translations,  also  as  described  in  section  3.1.  The 
third  parameter  of  the  image  function,  9  (see  figure  3.1)  is  an  encoding  of  the  orientation  of  the 
angle  passing  through  the  receptive  field.  We  discuss  how  this  encoding  is  developed  elsewhere  in 
this  document. 

Each  cell  A(p,  8),  A :  Rx  (0,  2-k]  — »  SR  in  the  accumulator  array  A  may  represent  a  line  or  an  edge 
in  the  image.  For  this  initial  work,  we  limit  the  exploration  to  detecting  lines.  The  variable  9  is 
the  same  9  as  mentioned  in  the  previous  paragraph,  and  parameterizes  the  slope  of  the  edge.  As 
in  Figure  4.1,  p  is  the  minimum  distance  of  this  edge  to  the  foveola2. 

The  accumulation  function  used  here  relies  on  a  critical  observation:  At  any  point  in  Cl,  both  the 
spatial  coordinates,  and  the  orientation  of  the  edge  are  known,  therefore  the  edge  is 
uniquely  determined.  For  this  reason,  only  one3  afferent  connection  is  required  from  a  Cl  cell  in 
VI,  Cl(x,y,9),  to  a  cell  in  the  accumulator,  A(p,9),  as  illustrated  in  Figure  4.1.  The  accumulator 
may  be  considered  as  a  two  dimensional  array  indexed  by  8  and  p.  It  is  straightforward  to  show 
that  for  any  point  x,  y  on  the  line, 


p  =  (xcos9  +  ys'm9)  (4.15) 

The  accumulator  cells  have  as  many  efferents  as  there  are  Cl  cells  in  the  corresponding  edge,  which 
could  be  on  the  order  of  a  thousand,  still  a  reasonable  biological  assumption. 

Here,  we  represent  a  straight  line  by  a  finite  set  of  points  (r,  9).  The  set  of  points  constituting  line  Z  is  defined  b^ 

A(p,9)  =  {r,9\  \p  —  t cos 9  —  r  sin  $  |  <d}  (4.16) 

where  5  is  some  “small”  positive  constant.  Thus,  the  accumulator  cell  with  parameters  (p,  (f>)  has 
a  value  equal  to 

n 

AP,e  =  J2C1*  (4-17) 


2  Although  the  foveola,  the  center  of  the  fovea,  has  a  non-infinitesimal  area,  it  is  quite  small,  and  the  term  is  used 
here  to  denote  the  center  of  the  foveola,  the  origin  of  the  coordinate  system 
2One  connection  is  not  robust,  however,  and  this  is  addressed  in  section  5.2 
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where  i  simply  denotes  a  cell  on  the  line  p,  6.  The  fact  that  the  orientation  returned  by  edge 
detectors  is  the  same  as  the  parameter  used  by  the  transform  potentially  provides  for  the  simple 
interconnection  architecture  mentioned  above. 
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Chapter  5 

Approach 


In  this  section,  we  describe  topics  addressed  in  the  initial  (STIR)  phase  of  the  project,  and  topics 
likely  to  be  addressed  in  the  second,  follow-on,  phase  of  the  project.  In  this  way,  the  reader  may 
see  the  entire  scope  of  the  effort. 

In  this,  the  “Approach”  section,  we  will  describe  experiments  run  and  also  the  results  of  those 
experiments. 


5.1  Estimating  Orientations 

Schoups  et  al.  [51]  indicate  that  there  is  no  particularly  preferred  orientation  in  the  visual  field 
of  monkeys.  Using  a  fixed,  small  number  of  Gabor-like  filters  suggests  a  reasonable  biological 
implementation,  just  a  few  orientation-sensitive  cells  at  every  point.  Yet,  we  must  ask,  if  we  only 
make,  say,  four  independent  measurements  of  orientation,  and  each  one  is  accurate  only  to,  say, 
one  part  in  thirty  (five  bits  precision),  how  can  eight-bit  precision  occur,  as  has  been  measured  in 
primates [51]?  In  Appendix  A,  we  show  how  four  measurements  of  the  orientation  of  a  particular 
edge  can  be  used  to  produce  an  accurate  estimate  of  the  actual  orientation.  However,  though 
that  appendix  demonstrates  that  a  linear  machine  can  determine  a  single  estimate  from  the  four 
measurements,  it  does  not  consider  in  detail  how  a  biological  collection  of  cells  could  perform  this 
operation. 

In  this  phase  of  the  research,  we  have  examined  the  method  described  in  the  appendix,  as  well  as 
a  number  of  other  methods  to  validate  how  well  they  work  on  long  and  short  line  segments  and 
how  tolerant  they  are  to  errors  in  the  original  measurements.  In  this  phase,  we  also  determine  the 
best  way  to  use  a  set  of  measurements  directly  to  select  the  cells  or  cell  neighborhood  in  A  to  be 
accumulated.  In  Figure  5.1,  we  illustrate  the  result  of  a  simulation  using  multiple  neurons  with 
only  a  few  bits  of  precision.  By  exploiting  consistency  (as  determined  by  the  accumulator),  high 
precision  results  can  be  obtained  from  these  low  precision  calculations. 

In  the  following  subsections,  we  explore  a  variety  of  options  for  estimating  the  orientation  of  an 
edge  from  a  limited  number  of  measurements.  We  assume  n  measurements  of  orientation  have  been 
made,  and  we  wish  to  find  a  best  estimate  of  the  actual  orientation.  The  interpolation  is  to  be 
accurate  over  a  range  of  0  to  it  which  has  been  quantized  into  m  possible  angles.  For  convenience, 
we  choose  m  to  be  180,  so  that  a  change  of  A 6  =  7r/180  corresponds  to  one  degree. 

Typically,  we  might  make  6  measurements  using  the  Gabor  filters,  and  interpolate  those  6  mea¬ 
surements  to  an  accuracy  of  one  degree. 

In  the  following  explanatory  sections,  we  will  use  one  of  the  simple  inputs  shown  below: 
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We  experimented  with  a  number  of  ways  in  increment  the  accumulator  based  on  a  set  of  directional 
measurements.  These  are  described  below: 

5.1.1  Estimating  Orientation  Using  the  Traditional  Hough  Transform  (mode  0) 

Mode  0  is  the  classical  Hough  Transform  (HT).  At  each  point,  the  entire  accumulator  is  incremented. 
First,  the  maximum  strength  of  the  response  is  tested,  and  if  it  is  significant,  a  plot  of  the  curve 
p  =  x  cos  8  +  ys'mO  is  generated  by  allowing  0  to  range  over  permissible  values.  We  implemented 
this  method  only  because  it  provides  a  reference  to  a  well-known  method.  Interestingly,  we  found 
the  Hough  Transform  to  not  only  be  slower,  but  less  accurate  than  other  methods,  probably  due 
to  the  fact  that  peaks  in  the  transform  space  tend  to  not  be  symmetric.  Experimentally,  it  works 
well,  but  other  estimators  are  faster,  and  the  HT  does  not  meet  the  requirement  to  be  biologically 
plausible.  Figure  5.3  illustrates  the  accumulator  produced  by  the  classical  HT,  and  the  image 
reconstructed. 

5.1.2  Estimating  Orientation  By  Interpolation  (mode  1) 

This  mode  uses  the  n  measurements  of  directional  derivative  magnitude  and  then  linearly  inter¬ 
polates  them  over  a  total  of  180  1-degree  measurements.  Figure  5.4  illustrates  the  results  on  the 
test  image.  Unfortunately,  this  method  seldom  works,  probably  due  to  the  idea  that  the  true  angle 
may  be  derived  from  a  single  linear  interpolation  is  simply  false.  This  method  may  be  compared 
with  mode  5,  which  uses  a  parabolic  interpolator  and  works  very  well. 

5.1.3  Estimating  Orientation  Using  the  Pseudo-inverse  (mode  2) 

In  Appendix  A,  we  illustrate  how  a  gradient  vector  may  be  estimated  from  4  noisy  estimates.  It  is 
straightforward  to  extend  this  to  a  larger  number  of  estimates  than  4,  and  we  did  so. 

A  single  noisy  gradient  vector,  projected  onto  all  four  directions  produces  four  (or  n  in  the  general 
case)  noisy  measurements.  The  method  then  returns  the  original  vector  quite  precisely.  However, 
if  the  n  measurements  are  not  simply  projections  (as  the  Gabor  does  not  return  strict  projections), 
the  estimated  value  may  be  quite  imprecise.  We  ended  up  unable  to  use  this  method  for  n  >  4. 


5.1.4  Estimating  Orientation  Using  the  maximum  Gabor  Output  (mode  3) 

At  each  point,  a  number,  n  of  filters  have  been  applied.  The  largest  response  is  chosen  as  the 
correct  angle.  The  problem  with  this  approach  is  illustrated  in  figure  5.2  which  shows  a  single  line 
which  has  been  sampled  by  a  collection  of  filters  several  times.  The  line  subtends  an  angle  of  30 
degrees  with  the  x  axis.  Assuming  there  are  four  Gabor  filters  (0,  45°,  90°,  135°).  In  each  case,  the 
45°  filter  will  have  the  strongest  output.  However,  if  no  interpolation  is  done,  those  four  responses 
will  not  align.  Figure  5.5  illustrates  the  nearly  perfect  results  on  the  test  image. 

5.1.5  Estimating  Orientation  by  Incrementing  all  the  Directions  (mode  4) 

Because  all  n  directions  are  incremented,  a  significant  amount  of  blur  allows  a  successful  detection. 
This  mode  is  attractive  because  no  interpolation  is  needed  at  all,  except  for  the  natural  sum¬ 
ming  of  the  measurements.  Furthermore,  this  is  the  mode  originally  envisioned  in  the  proposal. 
Surprisingly,  it  also  has  the  best  performance.  Figure  5.6  illustrates  the  results  on  the  test  image. 


5.1.6  Estimating  Orientation  by  Fitting  a  Parabola  (mode  5) 


First,  find  the  maximum  response  from  the  n.  Let  that  angle  be  denoted  62  and  the  filter  output 
at  that  point  be  1/2 ■  Let  the  angle  measurement  on  the  left  be  0\  with  output  y\  and  similarly  the 
point  on  the  right  be  63,1/3. 


Let  the  distance  (in  degrees)  be  7,  and 
passes  through  these  3  points  satisfies 

move  the  origin  to  allow  62  be  zero. 

A  parabola  which 

yi 

=  a(— -y)2  —  67  +  c 

(5.1) 

2/2 

=  c 

(5.2) 

2/3 

=  a(  7)2  +  67  +  c 

(5.3) 

This  system  of  equations  is  satisfied  by 

c 

=  2/2 

(5.4) 

b 

=  (2/3  -  2/1 ) / 27 

(5.5) 

a 

2/3  -  67  -  c 

72 

(5.6) 

Then,  the  location  of  the  strongest  response  is  found  by  differentiating  the  parabola  and  setting 

the  result  to  zero  producing 

9  =  —b/2a 

(5.7) 

This  approach  can  only  fail  in  the  degenerate  case  that  all  3  y’s  are  the  same.  In  that  case,  the 
estimated  direction  is  simply  that  of  62- 

Numerous  experiments  confirm  this  approach  works  best.  The  only  problem  is  the  division  of 
equation  5.7  is  challenging  to  motivate  in  a  biological  network.  Figure  5.7  illustrates  the  results  on 
the  test  image; 


5.1.7  Estimating  Orientation  by  Linearly  Interpolating  the  Strongest  Two  Di¬ 
rections  (mode  6) 

One  other  approach  is  to  choose  only  two  measurements,  the  one  having  the  strongest  response,  and 
one  of  its  neighbors,  and  performing  a  linear  interpolation.  This  approach  does  not  work  because 
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this  would  require  that  the  interpolation  be  less  than  the  already-measured  maximum.  Figure  5.8 
illustrates  the  results  on  the  test  image. 


5.2  Accumulation 

A  single  connection  from  the  image  plane  to  the  accumulator  raises  the  possibility  that  several 
nearby  points  in  the  accumulator  may  be  each  individually  incremented,  without  one  particular 
point  emerging  as  the  clear  representative.  In  the  traditional  Hough  transform,  this  is  most  fre¬ 
quently  dealt  with  by  blurring  the  accumulator. 

In  the  traditional  Hough  transform,  however,  the  accumulation  of  points  near  a  peak  becomes 
problematic  because  the  distribution  of  points  is  not  isotropic.  Figure  4.4  illustrates  a  traditional 
Hough  Transform  with  distinct  peaks.  Even  though  the  peaks  are  distinct,  near  the  peaks  the 
distribution  is  strongly  non-isotropic,  and  a  large-scale  smoothing  model  will  not  be  appropriate. 
We  attempted  to  solve  the  problem  of  finding  the  peak  by  inventing  a  new  clustering  algorithm 
which  we  called  BlackHole.  In  this  approach,  each  nonzero  point  in  the  accumulator  represents 
a  particle  in  free  space  with  a  mass  equal  to  its  brightness,  and  moving  in  response  to  the  net 
gravitational  field  of  all  the  points.  Gradually,  the  points  all  collapse  into  a  single  very  bright  spot 
which  will  be  the  cluster  center. 

Since  this  effort  was  a  bit  of  a  diversion  from  the  main  emphasis  of  this  project,  we  allocated 
only  one  person- week  to  it.  Unfortunately,  the  problem  turned  out  to  have  a  number  of  subtle 
details,  was  in  general  more  complex  than  we  had  anticipated,  and  after  one  person-week,  we  had 
to  abandon  the  effort  before  it  yielded  a  solution.  Fortunately,  modes  3-5  yield  tight,  dense 
clusters  which  are  relatively  round.  We  found  that  simple,  small  area  blurring,  followed  by  local 
maximum  detection  was  sufficient  for  modes  3,  4,  and  5. 
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Figure  5.1:  Graphs  show  average  precision  (measured  in  radians)  in  determining  the  orientation 
of  a  straight  line,  as  a  function  of  the  number  of  bits  of  accuracy  used  in  the  estimate.  From  top 
to  bottom,  the  curves  illustrate  the  effects  of  having  2,  3,  5,  and  10  points  along  that  line.  It 
is  clear  that  more  than  ten  points  do  not  increase  the  precision.  The  parameters  are  determined 
by  choosing  points  roughly  along  the  line  (“rough”  is  determined  by  the  bits  of  precision),  and 
then  measuring  the  p ,  <j>  parameters  of  the  line.  These  parameters  are  also  computed  only  to  the 
precision  of  the  bits  specified  along  the  horizontal  axis.  Since  the  possible  range  was  n r,  an  accuracy 
of  0.01  at  5  bits  is  equivalent  to  one  part  in  300,  or  over  eight  bits  of  precision. 


Figure  5.2:  If  the  orientation  resolution  is  too  small  (an  insufficient  number  of  directional  filters  is 
used),  the  sampled  points  will  not  be  collinear 
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Figure  5.3:  (MODE  0)  LEFT: Accumulator  resulting  from  finding  lines  in  the  test  image,  using  the 
traditional  Hough  Transform.  RIGHT:  The  image  reconstructed  from  the  accumulator.  Note  that 
the  output  image  (512  x  512)  is  larger  than  the  input  image  (256  x  256),  so  the  reconstructed  image 
shows  the  edges  in  the  upper  left  quadrant. 


Figure  5.4:  (MODE  1)  LEFT: Accumulator  resulting  from  finding  lines  in  the  test  image,  using 
mode  1,  which  does  not  work.  RIGHT:  The  image  reconstructed  from  the  accumulator. 
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Figure  5.5:  (MODE  3)  LEFT: Accumulator  resulting  from  finding  lines  in  the  test  image,  using 
mode  3,  which  works  well.  RIGHT:  The  image  reconstructed  from  the  accumulator.  Note  that  the 
output  image  (512  x  512)  is  larger  than  the  input  image  (256  x  256),  so  the  reconstructed  image 
shows  the  edges  in  the  upper  left  quadrant. 


Figure  5.6:  (MODE  4)  LEFT: Accumulator  resulting  from  finding  lines  in  the  test  image,  using 
mode  4,  which  works  well.  RIGHT:  The  image  reconstructed  from  the  accumulator.  Note  that  the 
output  image  (512  x  512)  is  larger  than  the  input  image  (256  x  256),  so  the  reconstructed  image 
shows  the  edges  in  the  upper  left  quadrant. 
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Figure  5.7:  (MODE  5)  LEFT: Accumulator  resulting  from  finding  lines  in  the  test  image,  using 
mode  5,  which  works  well.  RIGHT:  The  image  reconstructed  from  the  accumulator.  Note  that  the 
output  image  (512  x  512)  is  larger  than  the  input  image  (256  x  256),  so  the  reconstructed  image 
shows  the  edges  in  the  upper  left  quadrant. 


Figure  5.8:  (MODE  6)  LEFT: Accumulator  resulting  from  finding  lines  in  the  test  image,  using 
mode  6,  which  does  not  work  well.  RIGHT:  The  image  reconstructed  from  the  accumulator.  Note 
that  the  output  image  (512  x  512)  is  larger  than  the  input  image  (256  x  256),  so  the  reconstructed 
image  shows  the  edges  in  the  upper  left  quadrant. 


25 


Chapter  6 

Results 


Figure  6.1  illustrates  the  type  of  image  we  used  for  testing.  Instead  of  natural  images,  we  used  this 
class  of  test  images  because  it  is  possible  to  carefully  control  the  parameters  of  interest.  The  reader 
should  look  at  the  first  image  and  identify  the  most  distinctive  straight  line  in  that  image.  Most 
human  observers  choose  the  straight  line  at  an  angle  of  —45°  to  the  x-axis  in  the  lower  right.  That 
line  consists  of  five  microlines.  However,  there  are  also  two  other  lines  of  interest  in  this  image,  a 
horizontal  line  at  the  upper  left  ,  and  another  horizontal  line  about  1/3  of  the  way  down.  Although 
both  of  those  macrolines  are  only  of  length  2,  nonetheless,  they  are  collinear  segments,  and  we 
would  expect  a  program  which  models  human  behavior  to  detect  them. 

This  experiment  was  motivated  by  previous  work  in  neural  modeling  [67,  59,  30] ,  using  similar  test 
images.  Gintautas  et  al.  state  “we  employ  abstract  computer-generated  shapes  consisting  of  short, 
smooth  contour  segments  that  could  either  be  globally  aligned  to  form  wiggly,  nearly  closed  objects 
( amoebas ,  or  else  randomly  rotated  to  provide  a  background  of  locally  indistinguishable  contour 
fragments  ( clutter ).”  In  our  case,  we  seek  not  amoebas  but  lines,  but  the  philosophy  is  identical. 
By  using  test  images  like  this,  we  can  control  a  variety  of  parameters  such  as: 

length  Here,  the  term  length  refers  to  the  length  of  the  macroline.  It  is  the  total  number  of 
microlines  which  make  up  the  nracroline.  For  example,  the  length  of  the  45°  macroline 
referred  to  above  is  5. 

Microline  orientation  A  nracroline  consists  of  one  or  more  microlines,  and  each  microline  has 
an  orientation.  The  value  orientation  is  actually  the  variance  of  the  orientation  of  the  micro¬ 
lines  which  make  up  a  single  nracroline.  Performance  as  a  function  of  variance  of  nricroline 
orientation  is  discussed  in  section  6.5.3. 

We  have  tested  the  peak  detection  algorithm  extensively,  and  will  show  performance  results  later 
in  this  section.  First,  we  discuss  the  nature  of  the  data  and  the  experiments. 


6.1  The  Experiments 

Using  modes  3,  4,  and  5,  we  process  this  inrage,  find  the  accumulator,  and  reconstruct  the  ac¬ 
cumulator  to  produce  Figures  6.2  -  6.4.  The  brightness  of  the  line  in  the  reconstructed  inrage  is 
proportional  to  the  confidence  the  algorithm  has  in  the  detection. 

We  then  repeat  this  experiment  using  the  second  inrage.  Remember,  these  two  images  differ  only 
in  the  small  random  perturbation  applied  to  the  diagonal  line.  Interestingly,  all  three  methods 
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Figure  6.1:  Two  similar  images  from  the  set  of  test  images.  The  45°  line  in  the  lower  right  of  the 
left  image  is  deformed  by  a  5°  random  variation  to  form  the  image  on  the  right.  Note,  on  some 
computer  screens,  lines  may  disappear  due  to  the  coarse  sampling  required  to  show  both  images 
on  the  same  page  of  the  report.  Also  note  that  aliasing  effects  occur  in  these  images  which  were 
corrected  as  illustrated  in  Figure  8.1 


fail  to  identify  the  diagonal  line  as  having  high  likelihood.  We  have  run  experiments1  with  humans 
and  found  a  similar  result,  although  humans  will  sometimes  identify  the  dark  region  around  the 
diagonal  line  as  a  collinear  segment.  Even  small  variations  away  from  collinearity  caused  significant 
changes  in  the  human’s  sensitivity  to  lines. 


6.2  Detection  as  a  Function  of  Length 

Figure  6.8  illustrates  the  probability  of  correct  detection  as  a  function  of  length  of  the  nracroline. 
This  figure  is  an  average  over  a  variety  of  values  for  variation  in  nricroline  orientation  and  nracroline 
orientation. 


6.3  Detection  as  a  Function  of  Microline  Orientation 

Figure  6.9  shows  the  probability  as  a  function  of  variance  in  the  orientation  of  the  microlines.  It 
is  particularly  interesting  to  compare  the  performance  of  the  method  with  the  performance  of  a 
human  in  this  context.  For  example,  note  the  two  images  in  Figure  6.1  differ  only  by  the  random 
variation  of  the  microlines  in  the  single  45°  line,  and  that  variance  is  only  five  degrees.  Yet,  when 
the  simulation  is  run  on  the  left  image,  the  45°  line  is  perceived  as  the  most  distinctive,  but  when 
it  is  run  on  the  right  image,  the  horizontal  line  near  the  center  of  the  image  replaces  the  45°  line 
as  the  most  distinctive.  Several  humans  that  we  polled  agree  with  this  distinction,  suggesting  that 
the  model  is  consistent  with  human  behavior. 

1The  results  of  these  experiments  can  only  be  classified  as  anecdotal,  as  they  were  performed  without  an  expert 
in  human  cognitive  psychology  to  design  them. 


27 


Figure  6.2:  The  first  test  image  processed  with  mode  3. 
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Figure  6.3:  The  first  test  image  processed  with  mode  4. 
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Figure  6.4:  The  first  test  image  processed  with  mode  5. 
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Figure  6.6:  The  second  test  image  processed  with  mode  4. 
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probability  of  detection  based  on  length 


Figure  6.8:  Detection  probability  of  macrolines  as  a  function  of  length.  The  performance  is  lim¬ 
ited  by  random  variation  in  orientation  of  microlines.  See  discussion  in  text.  In  this  image,  all 
background  clutter  was  the  same  intensity  as  the  target  line. 


This  may  be  compared  with  the  work  of  Gintautas  et  al.,  [30]  who  work  with  “amoebas”  which 
are,  by  design,  figures  in  which  local  microlines  (our  terminology)  are  very  closely  aligned.  Lack  of 
local  alignment  should  produce  a  negative  decision. 

It  is  reasonably  easy  to  see  why  the  collinearity  requirement  is  so  sensitive.  Imagine  two  short, 
nonparallel  line  segments  one  with  orientation,  say  53  degrees  and  the  other  47  degrees.  Suppose 
their  centers  both  lie  on  the  same  50°  line.  Even  though  they  are  both  have  orientations  close 
to  that  of  the  line  on  which  they  lie,  when  one  extends  the  segments,  we  see  that  they  are  really 
quite  different  lines.  The  transform-based  detection  method  we  use  does  not  distinguish  distance 
between  the  line  segments,  so  they  are  not  considered  part  of  the  same  line.  It  is  clear  from 
other  experiments,  that  the  human  follows  a  curve,  if  the  segments  are  touching  or  very  close,  but 
apparently  not  if  the  segments  are  not  connected  as  in  this  case. 


6.4  Detection  as  a  Function  of  Clutter 

Not  surprisingly,  if  the  background  clutter  is  reduced,  the  probability  of  correct  detection  of  the 
true  nracroline  increases.  This  is  illustrated  in  Figure  6.10. 


6.5  Comparison  with  human  experiments 

In  this  section,  the  results  from  SLDSimulation  are  compared  with  results  from  [26]  and  [35]. 
Contour  detection  experiments  with  ‘yes’  and  ‘no’  responses  were  conducted  on  human  subjects  in 
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probability  of  detection  based  on  variation  of  micro  line  orientation 


Figure  6.9: 
microline. 


Detection  probability  of  macrolines  as  a  function  of  variance  in  the  orientation  of  the 
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probability  of  detection  based  on  length 


length 


Figure  6.10:  Detection  probability  as  a  function  of  length  for  a  low  signal-to-clutter  ratio  image. 
Here,  the  average  brightness  of  a  true  positive  pixel  is  double  that  of  a  clutter  pixel.  Compare 
with  Figure  6.8  in  which  the  clutter  and  signal  are  the  same  brightness.  This  is  consistent  with  the 
observation  at  the  beginning  of  the  report  that  even  very  subtle  lines  may  alert  the  line-detecting 
portion  of  the  brain  if  the  clutter  is  very  low. 
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[26]  and  [35]. 


6.5.1  Similarities  in  experimental  parameters 

Here,  we  use  the  terms  “microline”  and  “macroline”  to  distinguish  between  the  short,  sometimes 
collinear  segments  in  the  test  images  and  the  long  line  perceived  as  made  from  the  shorter  microlines. 
The  image  size  used  in  the  current  set  of  experiments  is  512  x  512,  and  the  image  is  divided  into 
squares  of  size  32  x  32.  Every  nricroline  was  constrained  to  fall  within  a  square.  The  straight  line 
was  constructed  first,  and  then  the  empty  squares  were  filled  up  with  microlines.  These  parameters 
are  similar  to  those  of  [26],  and  comparable  to  those  of  [35]. 

6.5.2  Differences  in  experimental  parameters 

Experiments  outlined  in  [26]  and  [35]  involve  detection  of  contours,  while  we  are  interested  in 
detection  of  straight  lines.  Line  segments  were  used  in  place  of  Gabor  elements.  The  position  of 
each  nricroline  was  constrained  to  fall  within  a  square,  and  the  position  within  the  square  varied 
in  Gaussian  manner  (with  high  probability  of  the  nricroline  falling  in  the  center  of  the  square), 
as  opposed  to  uniform  distribution  commonly  reported  in  literature.  Path  angle  remains  zero 
throughout,  as  we  are  interested  in  straight  lines.  It  varies  from  element  to  element  in  the  papers. 
Human  subjects  are  asked  “do  you  see  a  curve  (or  line)  at  all?”  We  ask  our  algorithm,  ”What  is 
the  orientation  of  the  dominant  line  you  see  (if  you  see  one)?” 


6.5.3  Results  Compared  with  Human  Performance  in  Literature 

Figure  6.11  shows  the  relation  between  variance  in  nricroline  orientation  along  a  curve  as  observed  by 
humans.  The  nricroline  orientation  is  varied  randomly  in  a  Gaussian  manner.  For  our  experiments, 
we  consider  only  straight  lines,  and  these  correspond  the  zero  angle  case,  (y  axis  of  figure).  It  can 
be  seen  that  for  the  zero  degree  case,  as  the  variation  increase  from  0  degrees  to  30  degrees,  the 
percentage  of  correct  detection  by  hunrans  decreases  from  100%  to  around  80%  at  15  degrees,  and 
to  almost  50%  at  30  degrees. The  same  trend  is  shown  by  SLD Simulation,  except  that  we  did  not 
ask  the  algorithm  “do  you  see  a  curve  at  all?” 


FIGURE  1 1.  Results  of  Expt  III.  The  solid  symbols  plot  the  results  for  orientations  of  the  element  with  respect  to  the  path 
(a)  of  ±  1 5  deg,  and  the  open  symbols  for  orientations  of  ±30 deg.  The  dashed  line  shows  the  results  from  Expt  I  (a  =  0). 


Figure  6.11:  Results  from  experiment  4  of  [26]. 

The  set  of  results  in  Figure  6.13  are  from  [35].  Positional  sigma  refers  to  the  sigma  of  the  2D  Gaus¬ 
sian  determining  the  position  of  the  microlines  along  the  path  of  the  straight  line,  but  constrained 
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probability  of  detection  based  on  variation  of  micro  line  orientation 


Figure  6.12:  Red:  results  from  Field  et  al.  showing  probability  of  detection  as  a  function  of  variance 
in  microline  orientation.  Blue:  SLDS  Simulation  results. 

to  lie  within  a  square.  We  are  primarily  interested  in  the  lower  left  figure  (unfilled  circles)  which 
reports  five  degree  path  angle  contours.  Hess  and  Dakin  [35]  do  not  report  zero  degree  behavior, 
so  5  degrees  is  as  close  as  available  for  comparison. 

Figure  6.14  compares  our  results  for  straight  lines  with  [35]  at  5  degrees. 

The  experiments  whose  results  are  in  Fig.  6.15  are  similar  to  those  in  figure  6.11  but  compare 
different  parameters. 

Orientation  sigma  refers  to  the  sigma  of  the  Gaussian  distributed  variation  of  microline  orientation 
from  the  straight  line  path.  The  difference  once  again  is  in  the  path  angle,  which  has  been  set  to 
zero  degrees  in  Fig.  6.16.  The  result  is  hence  most  close  to  foveal  detection  of  a  five  degree  path. 
Note  that  Figure  6.16  includes  not  only  our  results  in  blue  but  a  restatement  of  the  results  from 
Hess  [35].  If  the  length  of  the  straight  line  had  been  set  to  8  microlines,  as  is  the  case  in  [35],  we 
would  have  obtained  a  100%  detection  for  the  straight  line  by  SLDSimulation.  Reducing  the  line 
length  to  6  affected  the  detection  rate,  and  reduced  it  for  higher  sigma  in  the  orientation  variation 
value.  Further  reduction  in  the  line  length  would  lead  to  even  more  rapid  deterioration  in  detection 
rate. 
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Figure  6.13:  Positional  Sigma  results  from  [35]. 


probability  of  detection  based  on  variation  of  micro  line  position 


Figure  6.14:  Positional  Sigma  comparable  results. 
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Figure  6.15:  Orientation  Sigma  results  from  [35]. 


probability  of  detection  based  on  variation  of  micro  line  orientation 


Figure  6.16:  Orientation  Sigma  comparable  results. 
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Chapter  7 


Future  Work 


7.1  Correlation  with  Human  Behavior 

Our  experiments  in  comparing  human  perception  and  simulation  results  can  only  be  called  “anec¬ 
dotal,”  since  we  do  not  have  the  skill  or  facilities  to  design,  set  up,  and  execute  experiments  with 
humans.  We  need  to  present  a  human  with  a  set  of  experiments  similar  to  those  described  in  section 

6.1  and  compare  the  results. 

Further  work  on  this  project  will  require  human  experiments  with  an  expert  in  human  cognitive 
psychology.  Unfortunately,  there  do  not  seem  to  be  psychologists  near  NCSU  who  have  interest  or 
time;  e.g.  we  located  one  expert  at  Duke  University  who  has  the  skills  and  facilities,  but  he  has 
just  received  significant  miliary  funding  and  will  not  be  available  for  approximately  two  years.  He 
suggested  some  other  psychologists  in  other  states  who  have  the  requisite  skills.  However,  if  we 
cannot  identify  a  local  collaborator,  and  must  communicate  long-distance,  we  might  as  well  use 
people  we  know. 

We  have  had  a  collaboration  with  Dr.  David  Badcock  at  the  University  of  Western  Australia 
who  could  bring  to  the  project  an  expert-level  understanding  of  both  Computer  Vision  and  neu¬ 
rophysiology  of  vision,  while  practicing  as  a  cognitive  psychologist  with  a  huge  experience  base  in 
experimental  design.  Unfortunately,  travel  to  or  from  Australia  would  be  expensive.  It  might  be 
possible  to  perform  these  experiments  without  travel.  This  remains  to  be  investigated. 

7.2  Hardware 

The  proposed  architecture  requires  one  connection  from  a  given  Cl  cell  to  an  accumulator  cell.  Each 
accumulator  cell,  however,  may  receive  input  from  many  Cl  cells,  perhaps  thousands,  depending 
on  the  resolution [46].  We  propose  to  build  a  small  electronic  model.  First,  the  prototype  will  be 
designed  using  analog  hardware  functionally  equivalent  to  a  neural  network[14,  32].  The  design 
will  be  simulated  on  a  digital  computer  and  its  performance  evaluated.  Then,  if  feasible,  we  will 
construct  a  hardware  prototype  capable  of  processing  a  reasonably  (e.g.  512  x  512)  sized  image. 
Both  analog  and  FPGA  implementations  will  be  evaluated. 
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7.3  Learning 


The  architecture  proposed  above  assumes  an  interconnection  network  which  is  specific  to  detecting 
straight  lines  as  they  are  presented  in  the  visual  field.  One  would  suspect  that  such  hardware 
is  built-in,  since  detection  of  the  horizon  is  so  critical  to  self  orientation,  but  biological  evidence 
is  lacking  to  indicate  how  much  is  built-in  and  how  much  is  learned.  Here,  we  assume  the  basic 
architecture  is  built  in,  but  the  weights  must  be  learned.  This  philosophy  is  consistent  with  the  work 
of  Linsker  [43]  who  assumed  a  locally-connected  array  of  light  detectors  with  local  interconnections, 
and  showed  how  learning  based  on  the  Hebb  rule  can  produce  receptive  fields. 

The  Hebb  rule  defines  how  the  interaction  between  two  neurons  may  be  associated.  Given  two 
neurons,  say  x  and  y,  connected  so  that  x  provides  input  to  y,  the  synapse  between  the  neurons 
changes  according  to  At]  =  ryx/y.  where  x  denotes  the  firing  rate  of  neuron  x.  Thus,  if  firing 
neuron  x  causes  neuron  y  to  fire,  then,  when  x  fires  again,  y  is  (slightly)  more  likely  to  fire,  rj 
is  the  learning  rate.  Hebb  learning  has  been  described  in  many  different  variations  since  Hebb’s 
original  work  in  1949  [34] .  However,  when  compared  to  more  computer-appropriate  algorithms  such 
as  backpropagation,  it  remains  the  most  biologically-plausible.  Similarly,  we  hypothesize  that  Hebb 
learning  can  result  in  interconnection  weights  which  will  detect  long  straight  lines. 

In  this  component  of  the  research,  we  will  ask  how  it  is  possible  to  arrange  the  accumulator  array 
in  such  a  way  that  Hebb  learning  can  produce  the  interconnection  weights  which  in  turn  produce 
the  accumulation  function  and  the  peak  detection  function. 

We  propose  to  also  use  Spiking  Neural  Network(SNN)  models  for  simulation  of  the  connections 
between  complex  cells  and  the  accumulator  array.  SNNs  use  spike  trains  to  convey  information, 
as  opposed  to  neural  networks  which  use  single  values  [64].  The  learning  rule  used  is  called  spike 
timing  dependent  plasticity,  and  is  similar  to  Hebbian  learning.  SNNs  work  with  a  large  number 
of  neurons,  and  may  be  trained  with  algorithms  similar  to  backpropagation.  Field  programmable 
gate  arrays  may  be  used  for  their  hardware  implementation.  To  handle  the  large  number  of  afferent 
connections  from  Cl  cells,  a  time  division  multiplexing  architecture  may  be  used[31]. 

7.4  Impact  on  Models  for  the  Visual  Cortex 

Neurophysiological  experiments  have  functionally  described  much  of  the  mammalian  visual  system. 
The  lower  levels  are  known  well,  as  mentioned  earlier.  We  know,  for  example,  that  there  are  cells 
which  respond  consistently  to  specific  stimuli,  such  as  edges. 

Here,  we  should  observe  that  all  the  tests  described  above  use  line  detection  rather  than  rising  edge 
or  falling  edge  detectors. 

At  the  lower  levels,  we  also  know  that  the  cortex  is  roughly  retinotropic.  That  is,  cells  which  respond 
to  nearby  stimuli  in  the  visual  field  are  themselves  close  together  in  the  cortex.  However,  other  than 
those  rather  simple  relationships,  other  researchers  have  not  found  good  mapping  between  image 
features,  geometry,  and  computational  organization.  The  approach  described  here  identifies  very 
simple  features,  straight  lines,  but  it  potentially  provides  an  understanding  of  other  architectural 
features  in  the  vision  system. 

In  section  4.2  we  postulated  the  presence  of  a  collection  of  cells  which  function  as  accumulators  of 
evidence,  evidence  which  would  suggest  the  presence  of  one  or  more  straight  edges  in  the  scene. 
We  have  shown  that  the  interconnections  between  Cl  cells  and  accumulator  cells  is  not  impossible, 
but  we  have  not  demonstrated  biological  evidence  for  such  an  accumulator.  A  significant  objective 
of  the  proposed  research  is  to  address  the  question,  “If  such  an  accumulator  exists,  how  can  its 
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presence  be  identified  and  measured  in  a  biological  neural  network?”  Of  course,  we  do  not  currently 
have  an  answer  to  that  question,  but  a  few  directions  of  investigation  may  be  suggested  for  future 
work: 


•  Although  each  Cl  cell  is  only  required  to  have  one  afferent  to  the  accumulator,  each  accu¬ 
mulator  may  have  many  inputs.  Therefore,  to  keep  the  interconnection  density  down  in  the 
accumulator,  some  architectural  variations  may  be  provided  to  trade  off,  for  example,  res¬ 
olution  in  the  radial  direction  for  interconnection  density.  For  example,  if  the  accumulator 
were  represented  in  polar  coordinates,  the  position  of  image  features  would  be  measured  more 
accurately  for  features  close  to  the  center  of  the  visual  field.  Could  this  be  an  explanation  for 
the  well-known  “log-polar”  image  representation  which  has  been  detected  in  the  mammalian 
visual  system?  Of  course,  the  same  argument  for  varying  resolution  would  apply  to  pixels  as 
well  as  to  edges.  In  future  work,  it  might  be  possible  to  answer  this  question. 

•  Carl  Weirnan  published  a  number  of  papers  on  log-polar  transform  (See  section  3.2),  includ¬ 
ing  one[62]  which  is  particularly  relevant  to  this  proposal.  In  that  paper,  he  used  the  duality 
between  lines  and  points  to  show  that  if  one  simply  used  log  p,  9  as  the  parameters  of  the 
accumulator  instead  of  p,  9,  the  “log-Hough”  array  was  identical  to  the  log-polar  representa¬ 
tion  of  an  image  described  in  section  3.2.  We  propose  to  extend  this  understanding  to  take 
advantage  of  our  earlier  observation  that  using  orientation  knowledge,  one  can  reduce  the 
number  of  axons  between  the  image  and  the  accumulator.  Weirnan  suggests  that  both  may 
occupy  the  same  portion  of  the  cortex,  taking  care  to  maintain  separation  between  image 
data  and  accumulator  data.  We  speculate  that  this  may  lead  to  an  understanding  of  the 
principles  underlying  the  existence  of  the  log-polar  representation  in  the  cortex. 

•  Following  on  the  previous  topic,  we  question  how  such  an  accumulator  may  be  laid  out  in 
the  cortex.  We  always  draw  the  accumulator  as  if  it  were  physically  located  separate  from 
the  Cl  area,  but  that  is  simply  a  convenience  for  thinking  about  what  feeds  signals  to  what. 
There  is  no  particular  reason  for  such  a  separation,  and  it  is  likely  that  the  accumulator  cells 
are  interspersed  between  the  imaging  cells.  We  would  hope  to  investigate  both  architectures, 
and  evaluate  average  connection  length,  axon  density,  and  learning/ adaptation  mechanisms. 
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Chapter  8 

An  Algorithm  for  Drawing  Straight 
Lines 


To  perform  the  simulations  described  in  section  6.1,  it  is  necessary  to  draw  straight  lines.  As  is 
well-known,  it  is  difficult  to  draw  straight  lines  which  are  very  thin  (one  pixel)  without  distortions  . 
Figure  6.1  illustrates  these  phenomena,  as  the  human  detects  angled  lines  differently  from  horizontal 
or  vertical  lines.  It  turns  out  that  the  Gabor  filter  which  we  use  to  estimate  directional  derivatives 
also  responds  differently  to  diagonal  and  vertical/horizontal  lines.  Correcting  this  is  referred  to 
in  the  Computer  Graphics  community  as  “anti-aliasing.”  The  paper  by  Xiaolin  Wu  [65]  probably 
represents  the  state-of-the-art  in  making  such  images  look  “good”  to  human  observers. 

However,  in  our  work  with  Gabor  filters,  we  found  a  sensitivity  that  humans  do  not  seem  to  have. 
Consider  a  white  horizontal  line  of  length  n  pixels  on  a  black  background,  and  think  of  it  as  a 
rectangle  1  pixel  high  and  n  pixels  long.  We  compare  that  with  a  diagonal  (45°)  line,  also  n  pixels 
long.  Place  a  circular  receptive  field  of  diameter  m  pixels  (m  >  2)  about  both  of  these  lines.  Within 
the  receptive  field,  there  is  more  “whiteness”  within  the  receptive  field  for  the  horizontal  line  than 
for  the  diagonal  line,  thus  the  Gabor  filter  will  produce  a  stronger  response. 

Our  solution  to  this  version  of  anti-aliasing  is  as  follows.  To  draw  a  line  of  length  1,  with  center  at 
point  xo,yo  and  orientation  9, 

•  Imagine  each  pixel  as  made  up  of  a  9  x  9  array  of  points,  equally  spaced  within  the  lxl  area 
of  the  pixel.  Refer  to  these  points  as  “subpixels.” 

•  Arrange  the  pixels  vertically  along  the  x  =  0  axis. 

•  For  each  subpixel,  rotate  it  by  6  about  the  center  of  the  line,  and  translate  the  resulting 
coordinates  by  xo,2/o- 

•  Determine  into  which  pixel  the  rotated-and-translated  subpixel  lies,  and  increment  the  bright¬ 
ness  of  that  pixel  by  1/81. 

This  correction  improves  significantly  the  uniformity  of  Gabor  for  lines.  Figure  8.1  below,  created 
using  the  new  method,  should  be  compared  with  Figure  6.1  which  creates  straight  lines  using  simple 
dot-plotting. 
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Figure  8.1:  A  collection  of  random  lines  segmented  created  using  the  new  anti-aliasing  algorithm. 
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Chapter  9 

Conclusion 


We  have  described  an  architecture  for  detecting  linear  features  in  images  which  is  biologically- 
plausible.  We  have  demonstrated  reasonable  performance  estimates  for  this  architecture  through 
simulation.  Simulation  has  confirmed  that  high  resolution  accumulation  can  be  implemented  with 
only  a  few  directional  derivatives,  implemented  using  Gabor  functions,  and  is  therefore  a  feasible 
model  for  how  humans  detect  straight  lines. 

If  the  simple,  sum-of-filter  output  (our  mode  4)  is  used,  the  accumulator  peaks  are  roughly  circularly 
symmetric,  and  have  easy-to-identify  peaks. 

Our  experiments  used  collinear  or  nearly-colinear  collections  of  short,  disjoint  line  segments.  Dis¬ 
joint  line  segments  are  much  more  strongly  identified  as  belonging  to  a  single  longer  line  when  the 
segments  are  very  close  to  collinear,  within  two  degrees,  and  the  computer  simulation  has  provided 
an  explanation  for  this  phenomenon  in  humans. 
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Appendix  A 

Determining  Orientation  from  Four 
Estimates 


At  the  SI  level  of  the  standard  model,  we  assume  we  have  four  estimates  of  the  orientation  of  an 
edge,  determined  using  four  different  Gabor  filters.  Here  we  assume  that  the  angle  of  observation 
is  or  can  be  rotated  into  a  frame  in  which  two  of  those  filters  are  orthogonal  and  we  denote  those  as 
vertical  and  horizontal ,  denoting  them  as  basis  vectors  [0  1]T  and  [1  0]r  respectively.  The  remaining 
two  filters  are  also  mutually  orthogonal,  and  are  represented  (in  the  frame  of  the  other  two  -  the 
base  frame)  by  [  (71  ]T  and  [  (7f  —  (71  ]T  respectively.  The  direction  of  the  edge  in  the  base 

frame  is  the  unknown  ordered  pair  G  =  [Gx  Gy]T ,  and  this  direction  has  been  measured  four  times, 
producing  /i  =  [  m  H2  fMi  IM  }T- 

Constructing  the  matrix  B  from  the  four  basis  vectors  produces  the  relation 

BG  =  n  ,  (A.l) 


Using  the  pseudoinverse  to  find  the  minimum  squared  estimate  of  G  for  this  overdetermined  prob¬ 
lem, 

G  =  {BT  B)~l  BT  ii  (A. 2) 

Once  we  have  the  vector  G,  those  two  numbers  define  the  direction. 


(BtB)~1Bt 


0  l 

3  0 


1 

2^2 

1(71 


i 

2U2 

_  1(71 


(A.3) 


Thus,  from  four  measurements  of  edge  direction,  a  least-squares  estimate  can  be  derived  by  a  simple 
sum-of-products . 

The  two  values  of  G  can  then  be  fed  into  a  quantization  neural  network  to  obtain  a  single  signal 
which  indicates  the  direction  (angle) .  The  quantization  neural  network  is  basically  a  straightforward 
lookup  table,  with  n  inputs  and  2"  outputs,  selecting  one  output  axon  as  a  function  of  the  code  on 
the  inputs;  functionally  a  binary  decoder. 
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Appendix  B 

The  Log-Polar  Transform 


The  material  in  this  appendix  was  written  by  Theju  Jacob  as  a  supplementary  projected  related 
to  this  project. 

B.l  Introduction 

In  this  report  we  examine  the  nature  of  computation  in  the  visual  cortex  as  discussed  in  literature. 
A  log  polar  transform  is  believed  to  be  the  mathematical  mapping  which  can  explain  how  signals 
received  by  the  eye  are  interpreted  by  the  human  brain.  We  start  with  a  brief  introduction  to 
various  elements  in  the  human  visual  pathway.  Beginning  with  experiments  which  provided  the 
first  clues  to  such  a  mapping,  we  trace  literature  which  first  put  forth  the  log  polar  idea,  and  then  at 
various  modifications  made  to  the  original  idea  when  more  experimental  results  became  available. 
We  also  look  at  the  two  major  explanations  given  for  such  a  structure  in  the  brain,  and  at  the 
arguments  for  and  against  each.  Section  2  looks  at  literature  primarily  from  the  1960s  and  1970s. 
Papers  from  1980s,  1990s  and  2000s  are  examined  in  section  3.  We  ponder  the  question  of  which 
of  the  two  explanations  might  be  more  accurate,  and  make  our  concluding  statements  in  section  4. 
The  human  visual  system  enables  us  to  process  and  respond  to  all  kinds  of  visual  stimuli.  Various 
structures  along  the  visual  pathway,  shown  in  Fig.l,  constitute  the  visual  system.  Starting  at  the 
retina,  neurotransmitters  travel  through  the  lateral  geniculate  nucleus  to  reach  the  visual  cortex. 
The  retina  itself  is  a  complex  structure,  with  rods  and  cones  acting  as  photodetectors,  and  layers  of 
horizontal,  bipolar,  amacrine  and  ganglion  cells  above  them  before  optic  nerve  fibres  are  reached, 
as  in  Fig. 2.  We  have  the  fovea  located  almost  at  the  center  of  the  retina,  representing  about  5°  of 
visual  angle.  The  highest  concentration  of  cone  cells  and  a  near  absence  of  rod  cells  are  found  at 
the  fovea.  The  outer  region  surrounding  the  fovea  is  called  perifovea. 

The  visual  cortex  is  divided  into  several  areas  named  VI,  V2,  V3  and  V4,  Fig. 3.  In  this  study, 
we  are  primarly  concerned  with  VI,  also  called  the  primary  visual  cortex  or  striate  cortex,  though 
other  areas  do  make  an  appearance  at  various  instances. 

B.2  Beginnings  -  1960s, 1970s 

The  widely  accepted  notion  today  is  that  the  projection  from  the  mammalian  retina  to  the  visual 
cortex  takes  the  form  of  a  logarithmic  conformal  mapping.  A  conformal  mapping  is  a  transformation 
that  preserves  local  angles.  The  idea  of  a  logarithmic  mapping  appeared  in  literature  for  the  first 
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Visual  cortex 


Figure  B.l:  Visual  Pathway  Structure  from  [28] 


time  in  1977[52] .  The  results  of  physiological  experiments  from  previous  decades  are  outlined  and 
the  reasons  why  a  logarithmic  mapping  explains  the  results  are  put  forward  by  Schwartz  in  [52]. 
Results  from  [18],  [39],  [40],  [1]  and  [2]  are  mentioned  extensively  in  [52].  In  [18],  Daniel  and  Whit- 
teridge  outline  experiments  involving  insertion  of  electrodes  into  brains  of  monkeys  and  baboons, 
and  discuss  the  responses  to  various  visual  stimuli.  The  authors  plotted  the  paths  traced  in  the 
visual  field  by  the  electrodes  in  various  cortex  locations.  The  magnification  factor,  which  is  the 
ratio  of  linear  distance  in  mm  between  two  points  on  the  cortex  to  the  angular  distance  in  degrees 
of  corresponding  points  in  the  visual  field,  was  computed  and  plotted  for  different  segments  of  the 
visual  field.  It  was  found  that  the  magnification  factor  was  the  same  for  all  points  with  the  same 
radius  in  the  visual  field,  regardless  of  the  angular  coordinates.  A  plot  of  the  results  obtained  by 
the  authors  is  shown  in  Fig. 4,  where  eccentricity  is  the  angular  distance  in  degrees  from  the  center 
of  the  visual  field.  We  can  see  that  the  magnification  factor  decreases  as  eccentricity  increases. 
However,  the  relationship  is  clearly  not  linear. 

In  [39]  and  [40],  further  studies  of  the  monkey  striate  cortex  are  outlined  by  Hubei  and  Wiesel. 
They  detail  orientation  columns  and  ocular  dominance  columns  and  their  mutual  relationship. 
Orientation  columns  refer  to  groups  of  cells  located  close  together,  which  show  a  preference  for 
a  particular  orientation,  while  ocular  dominance  columns  refer  to  groups  of  cells  which  show  a 
preference  for  a  particular  eye.  The  authors  inserted  probes  into  the  cortex  and  noted  how  the 
orientation  and  ocular  dominance  varied  across  the  cortex.  They  note  that  the  total  width  between 
an  array  of  orientation  slabs  that  take  in  all  of  180  degrees  of  orientation  and  a  left  plus  right  set  of 
ocular  dominance  columns  amount  to  0.5  to  1  mm  in  the  cortex.  Various  plots  outlining  the  changes 
are  outlined  in  their  work.  A  plot  relevant  to  our  current  investigation,  showing  the  variation  in 


49 


Cones 


Rods 


Figure  B.2:  Cells  in  Retina  from  [28] 


the  receptive  field  size  and  magnification  factor  with  distance  from  the  fovea  is  shown  in  Fig. 5.  We 
note  that  the  receptive  field  size  as  well  as  the  inverse  of  the  magnification  factor  increase  with  an 
increase  in  the  distance  from  the  fovea.  A  more  detailed  explanation  of  these  results  and  other 
results  related  to  the  mammalian  visual  system  can  be  found  online  in  [37]. 

In  [1]  and  [2],  Allman  and  Kaas  look  at  the  representation  of  the  visual  field  in  V2  as  well  as  in  the 
medial  area  of  the  visual  cortex  in  owl  monkeys.  Recordings  were  done  inserting  microelectrodes 
into  the  owl  monkey  cortex  and  stimulating  the  eye  using  bars  of  light.  They  discovered  that  the 
central  7°  of  the  horizontal  meridian  correspond  to  a  zone  of  the  cortex  across  the  width  of  V2, 
after  which  the  horizontal  meridian  representation  splits  and  forms  the  anterior  border  of  V2.  The 
vertical  meridian  was  found  to  be  covered  by  the  posterior  border  of  V2.  They  also  found  the  areas 
in  V2  which  represented  the  central  sections  of  the  visual  quadrant.  The  medial  area  was  found  to 
devote  a  greater  proportion  of  its  area  to  the  peripheral  areas  of  the  visual  field. 

The  results  from  Daniel  and  Whitteridge[18]  were  further  analyzed  by  Schwartz  [52],  and  the 
following  mathematical  expression  was  given  for  cortical  magnification: 

nr  =  k/r 

where  m  is  the  magnification,  k  is  a  constant,  and  r  represents  eccentricity  from  foveal  representation.  [52] 
states  further  that  magnification  is  a  differential  quantity,  and  that  we  are  interested  in  an  analytic 
function  whose  derivative  is  radially  symmetric  and  is  proportional  to  1/r.  The  following  complex 
logarithmic  function  is  put  forth  as  satisfying  the  before  stated  requirements. 

w  =  log  (re1^) 
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Figure  B.3:  Visual  areas  from  [44] 


where  (r,</>)  represents  a  point  in  the  visual  plane,  and  complex  number  w  represents  a  point  in  the 
cortical  plane. 

The  magnification  factor  for  the  lateral  geniculate  nucleus  or  LGN  takes  the  same  form  as  that  for 
the  striate  cortex.  The  density  of  retinal  ganglion  cells  as  well  as  reorganization  of  neurons  while 
projecting  from  LGN  to  striate  cortex  are  both  thought  to  bring  about  the  final  form  of  the  striate 
cortex  retinotopic  map.  The  structure  in  the  secondary  visual  cortex  is  also  shown  to  be  described 
by  a  logarithmic  conformal  mapping.  The  authors  also  discuss  how  such  a  mapping  can  lead  to  the 
structural  regularity  observed  in  the  striate  cortex  by  previous  works  [39]  and  [40].  In  particular, 
equal  angular  steps  in  the  visual  plane  must  transform  to  equal  linear  steps  in  the  cortex,  and  the 
complex  logarithm  is  shown  to  accomplish  such  a  transformation. 

In  [63],  log  spiral  grids  are  put  forth  for  digitization  of  images.  Rotation  and  magnification  of 
images  become  simple  pixel  shifts  on  using  the  described  technique,  thus  saving  computation  steps 
and  memory.  The  authors  describe  various  methods  to  achieve  the  necessary  transformation  of  the 
image  grids  by  introducing  various  tessellations. 

B.3  1980s, 1990s, 2000s 

In  this  section,  we  continue  with  the  investigation  of  the  presence  of  log  polar  mapping  in  the  human 
visual  cortex.  In  [53]  Schwartz  states  that  the  local  columnar  structure  as  well  as  the  retinotopic 
mapping  found  in  the  cortex  can  be  explained  by  means  of  complex  logarithmic  mapping.  The 
author  also  suggests  a  logarithmic  mapping  as  a  mechanism  for  size  and  rotation  invariance  in 
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Text-fig.  4.  The  magnification  factor  plotted  against  eccentricity.  The  data  have 
been  grouped  in  sectors  each  of  30°. 

Figure  B.4:  Plot  of  a  set  of  results  from  [18] 

vision.  The  expression  derived  previously  in  literature  has  now  been  modified  to  compensate  for 
the  divergence  at  zero  to: 

w  =  log(re*^  +  a) 

where  a  is  a  constant.  The  modified  expression  gives  us  a  linear  map  for  small  values  of  re ^  (< 
a)  and  a  logarithmic  map  for  larger  values.  The  author  also  demonstrates  how  changes  in  size  or 
rotation  of  objects  reduce  to  a  linear  shift  under  the  mapping. 

Schwartz  et  al.  in  [54]  studied  the  role  of  inferior  temporal  cortex  rather  than  the  primary  visual 
cortex.  The  IT  cortex  is  thought  to  play  a  crucial  role  in  visual  object  recognition,  by  enabling 
us  to  identify  shapes  and  patterns.  IT  neurons  have  large  receptive  fields  that  extend  to  both 
half  visual  fields,  and  includes  the  center  of  gaze.  Five  macaque  monkeys  were  tested  with  slits  of 
light,  complex  objects  and  Fourier  Descriptor (FD)  [3], [66]  stimuli.  The  FD  for  a  particular  shape 
was  obtained  by  expanding  the  boundary  orientation  function  as  a  Fourier  series.  The  boundary 
orientation  function  of  a  shape  in  turn  is  determined  by  measuring  the  orientation  of  the  shape’s 
boundary  at  regular  intervals  around  the  perimeter.  Using  FD  to  describe  a  shape  would  make  it 
invariant  to  position  and  shape  of  the  stimulus.  For  the  majority  of  IT  neurons  thus  tested,  size, 
position  and  contrast  variations  did  not  affect  the  output,  which  appeared  to  suggest  that  the  cells 
must  be  sensitive  to  the  overall  shape  of  the  stimulus.  The  authors  also  suggest  that  coding  of 
boundary  curvature  may  not  be  the  exclusive  function  of  IT  neurons,  as  some  of  them  have  been 
found  to  be  selective  for  color,  texture  and  spatial  frequency. 

The  striate  cortex  responses  to  visual  stimuli  in  rhesus  monkeys  are  recorded  in  [22]  by  Dow  et 
al.,  and  the  relationships  between  eccentricity  and  magnification  factor,  receptive  field  size  and 
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DISTANCE  FROM  FOVEA  (°) 

Fig.  6A  Graph  of  average  field  size  (crosses)  and  magnification  3  (open  circles)  against 
eccentricity,  for  five  cortical  locations.  Points  for  4°,  8°,  18°  and  22°  were  from  one  monkey; 
for  1  °,  from  a  second.  Field  size  was  determined  by  averaging  the  fields  at  each  eccentricity, 
estimating  size  from  (length  X  width)0-5. 

Figure  B.5:  Plot  of  a  set  of  results  from  [40] 


the  minimum  angle  of  resolution  are  examined.  The  authors  found  that  their  results  agreed  with 
previous  literature  [18], [39], [40]  with  minor  modifications,  except  that  they  did  not  find  inverse 
magnification  to  be  related  to  the  minimum  angle  of  resolution  by  a  constant  of  proportionality, 
like  in  [18].  They  also  note  that  the  parallel  relationship  between  field  size  and  inverse  magnifica¬ 
tion  do  not  hold  well  in  the  fovea,  as  inverse  magnification  becomes  smaller  than  field  size  at  an 
eccentricity  less  than  5°.  The  same  authors  present  a  detailed  mapping  of  the  monkey  striate  cortex 
in  [23] .  They  describe  a  new  technique  for  transferring  various  penetration  sites  on  the  cortex  onto 
a  map,  so  as  to  form  a  coherent  picture  of  the  cortex  and  its  various  areas.  Their  technique  enabled 
detailed  comparisons  of  different  hemispheres  of  the  cortex.  Parameters  measured  include  verti¬ 
cal/horizontal  magnification  anisotropy,  which  was  found  to  be  1.5:1  at  central  eccentricities,  and 
foveal  magnification,  which  was  found  to  be  different  along  different  meridians  in  different  fields. 

B.3.1  Emphasis  of  Central  Vision  in  Retina 

The  hypothesis  that  higher  retinal  ganglion  cell  densities  lead  to  the  final  form  of  striate  cortex 
retinotopic  map  was  previously  seen  in  literature  [52].  A  strong  argument  for  higher  density  of 
ganglion  cells  in  the  fovea  leading  to  higher  emphasis  on  central  vision  in  the  cortex  is  put  forth 
in  [61]  by  Wassle  et  al.  Those  authors  used  improved  techniques  for  determining  ganglion  cell 
densities  in  monkeys.  Amacrine  cells,  known  to  be  displaced  regularly  into  the  ganglion  cell  layer, 
were  discounted.  Results  are  shown  in  Fig. 6.  The  density  of  cones  were  found  to  decrease  from 
a  peak  of  250,000/mm2  at  the  fovea  to  11,500/mm2  at  3  mm  eccentricity.  Cone  pedicles,  the 
synaptic  terminals  of  cones,  are  nearly  absent  at  the  fovea,  but  increases  sharply  and  flattens  at 
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Figure  B.6:  Plot  of  a  set  of  results  from  [61].  From  [61]  -  ‘Density  gradients  of  cones,  cone  pedicles 
and  ganglion  cells  along  the  temporal  horizontal  meridian  of  a  macaque  monkey  fovea.  The  inset 
in  (A)  shows  the  retinal  area  examined.  The  blind  spot  (open  circle),  the  fovea(asterisk),  and  some 
blood  vessels  characterize  the  central  retina.  A  rectangl  3.5  x  2.5  mm  wide  was  cut  from  the  eye 
cup,  vertical  sections  were  cut  through  the  upper  part,  and  horizontal  sections  through  the  lower 
part.  (A)Cone  density;  (B)  cone  pedicle  density;  (C)  ganglion  cell  density;  (D)  comparison  of  the 
density  gradients  (modified  from  Wassle  et  al.,  1989).’ 


500/xrn,  peaking  at  a  plateau  of  32,000/mm2.  Ganglion  cells  have  a  density  profile  similar  to  that 
of  cone  pedicles,  but  more  ganglion  cells  are  present  at  lower  eccentricities.  A  plateau  is  present  at 
800/rrn,  with  a  maximal  density  of  60,000/mm2.  The  authors  point  out  that  the  change  in  ganglion 
cell  density  follows  a  profile  similar  to  change  in  magnification  between  central  and  peripheral  visual 
field  in  the  cortex,  and  hence  they  claim  that  ganglion  cell  density  can  fully  explain  the  variation 
in  the  cortical  magnification  factor. 

B.3.2  Emphasis  of  Central  Vision  in  Retino  Cortical  Pathway 

Van  Essen  et  al.[60]  undertook  a  detailed  study  of  the  visual  cortex  in  macaque  monkeys.  A  smooth 
variation  of  recording  sites  within  the  cortex  was  found  to  correspond  to  a  smooth  progression  in 
receptive  fields,  thus  preserving  all  neighborhood  relationships.  This  further  implied  that  the  cortex 
is  retinotopic.  Large  variations  were  found  in  the  relative  emphasis  on  different  parts  of  the  visual 
field  representation  in  striate  cortex  between  individuals.  The  authors  also  state  that  the  striate 
cortex  lacks  radial  symmetry  in  2  respects  -  greater  emphasis  is  placed  on  the  horizontal  meridian 
when  compared  to  the  vertical  meridian,  and  the  inferior  parts  of  the  visual  field  are  slightly  more 
emphasized  than  the  superior  parts  of  the  visual  field.  Plots  for  magnification  factors  vs  eccentricity 
were  obtained,  as  shown  in  Fig. 7.  The  authors  also  gave  new  expressions  for  linear  and  areal  cortical 
magnification  factors  which  better  fit  the  data  than  those  previously  available  in  literature.  For 
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example,  the  areal  cortical  magnification  factor,  the  ratio  of  area  in  the  cortex  in  mm2  to  the 
corresponding  area  in  the  visual  field  in  deg2,  was  now  given  by 

Ma  =  140(0.78  +  E)“2'20mm2/deg2 

where  E  is  the  eccentricity.  The  equation  for  Ma  predicted  a  ratio  of  4400  in  areal  magnification  at 
1°  vs  80°  eccentricity.  On  the  other  hand,  retinal  ganglion  cell  density  at  1°  was  found  to  be  about 
50  times  that  of  the  density  at  80°,  which  implied  that  central  vision  is  emphasized  much  more  in 
the  cortex  than  in  the  retina. 


c 


Figure  B.7:  Plot  of  a  set  of  results  from  [60].  From  [60]  -  ‘Cortical  magnification  as  a  function 
of  eccentricity.  Linear  magnification  (mm/deg)  along  iso-polar  contoursis  shown  in  (A)  and  along 
iso-eccentricity  contours  in  (B).  Areal  magnification  is  shown  in  (C).  Calculation  of  magnification 
factors  was  based  on  the  topographic  contours  illustrated  in  Fig.  5.  The  map  coordinates  of  the 
intersections  of  iso-eccentricity  and  iso-polar  contours  were  entered  into  the  computer  via  a  graphics 
tablet.  In  regions  where  there  was  sufficiently  detailed  information  about  visual  topography,  the 
spacing  of  contours  was  twice  as  close  as  that  shown  in  the  preceeding  figure.  Regression  lines 
were  calculated  according  to  the  least-squares  method.  For  each  compartment  the  length  of  each 
side  was  divided  by  the  length  of  the  corresponding  visual  field  compartment  to  obtain  the  linear 
magnification  factors  (Me  and  Mp  ,  and  the  area  of  the  cortical  compartment  was  divided  by  the 
area  of  the  visual  field  compartment  to  obtain  the  areal  magnification  factor.’ 


Azzopardi  and  Cowey[4]  further  supported  the  idea  that  the  central  visual  field  is  emphasized 
more  in  the  cortex  than  in  the  retina.  The  authors  used  a  retrograde  transneural  tracer  on  rhesus 
monkeys,  and  the  number  of  ganglion  cells  projecting  to  marked  areas  of  the  cortex  were  counted. 
They  computed  the  average  area  in  the  visual  cortex  receiving  a  projection  from  a  single  ganglion 
cell,  and  used  that  as  an  index  of  cortical  allocation.  If  cortical  allocation  were  to  mirror  the 
ganglion  cell  distribution  across  the  retina,  the  ratio  of  perifoveal  to  peripheral  cortical  allocation 
should  be  1.  However,  the  ratios  were  found  to  lie  in  the  range  from  3.29  to  5.93.  Therefore,  the 
authors  state,  the  emphasis  on  central  vision  in  the  cortex  cannot  be  explained  by  the  changes  in 
retinal  ganglion  cell  density  alone.  They  hence  contradict  the  results  from  [61]. 

The  same  authors  conducted  more  detailed  experiments  on  macaque  monkeys  and  presented  their 
results  in  [5].  This  time,  they  studied  the  allocation  of  neurons  for  the  fovea  in  the  dorsal  lateral 
geniculate  nucleus  as  well.  Once  again,  cortical  allocation  values  larger  than  1  were  obtained, 
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thereby  implying  that  the  cortical  map  of  the  retina  is  not  peripherally  scaled,  but  that  relatively 
more  cortex  per  ganglion  cell  is  alloted  in  the  fovea  and  surrounding  retina.  [5]  also  states  why 
arguments  about  estimates  of  ganglion  cells  representing  the  fovea  being  erroneous  do  not  hold  in 
the  current  set  of  experiments  -  the  authors  are  comparing  cortical  areas  with  their  total  number 
of  ganglion  cells,  including  laterally  displaced  ganglion  cells.  Arguments  for  higher  ganglion  cell 
counts  near  the  fovea,  as  well  as  for  the  presence  of  displaced  amacrine  cells  in  the  ganglion  cell  layer 
thereby  leading  to  erroneous  ganglion  cell  count  have  also  been  discounted.  The  results  obtained 
by  the  authors  indicate  that  the  representation  of  the  fovea  is  magnified  when  moving  from  the 
retina  to  the  lateral  geniculate  nucleus  and  from  the  lateral  geniculate  nucleus  to  the  striate  cortex. 
Further  support  to  the  idea  of  an  overrepresentation  of  the  fovea  in  the  cortex  is  given  by  Popovic  et 
al[48] .  The  authors  analyzed  data  from  already  available  literature  and  computed  the  relationships 
between  cortical  magnification  factor,  effective  ganglion  cell  separation  and  the  minimum  angle  of 
resolution.  A  set  of  results  are  shown  in  Fig. 8.  A  linear  relationship  was  found  between  effective 
ganglion  cell  separation  and  minimum  angle  of  resolution.  They  show  that  the  inverse  magnifica¬ 
tion  factor  and  the  before-mentioned  quantities  of  ganglion  cell  separation  and  minimum  angle  of 
resolution  are  also  related  linearly.  The  results  point  out  that  a  markedly  larger  amount  of  striate 
cortex  per  cell  and  amount  of  visual  cortex  is  devoted  to  fovea  and  surrounding  retina.  The  paper 
concludes  that  the  overrepresentation  of  the  fovea  and  immediately  surrounding  retina  are  caused 
by  an  additional  magnification  in  retino  cortical  pathway. 


Figure  B.8:  Plot  of  a  set  of  results  from  [48].  From  [48]  -  ‘Plot  of  the  product  of  effective  GC 
separation  (S,deg)  and  the  linear  cortical  magnification  factor  (M, mm/deg)  vs.  eccentricity  (filled 
symbols),  as  well  as  the  product  of  resolution  thresholds  (MAR,  deg)  and  M  vs.  eccentricity  (open 
symbols),  for  the  data  of  Engel  et  al.(1997)  Sereno  et  al.(1995),  describing  the  reserved  amount  of 
visual  cortex  (in  mm)  per  GC  and  the  amount  of  visual  cortex  needed  to  process  a  given  resolution 
threshold,  respectively.’ 


B.3.3  Later  studies,  use  of  fMRI 

The  question  of  how  information  transfer  happens  between  neurons  and  how  it  is  influenced  by 
cortical  organization  is  examined  in  [58].  There  are  2  possibilities  for  how  information  transfer 
occurs  between  neurons:  information  is  present  in  the  spike  rate  of  neurons,  or  in  the  precise 
timing  of  individual  spikes.  The  authors  examined  excitatory  post  synaptic  potentials  (EPSPs)  and 
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inhibitory  post  synaptic  potentials  (IPSPs)  of  neurons  and  the  question  of  how  they  are  balanced. 
For  orientation  selective  neurons,  for  example,  EPSPs  and  IPSPs  occur  frequently  for  optimally 
oriented  stimuli.  As  excitatory  inputs  bombard  a  cell  in  the  hundreds,  much  more  than  necessary  to 
depolarize  a  cell  membrane  from  resting  potential  to  spike  threshold,  the  authors  propose  that  the 
primary  role  of  inhibition  is  to  prevent  the  cell  firing  rate  from  reaching  saturation.  As  excitatory 
synapses  outnumber  inhibitory  synapses  by  6:1,  the  authors  suggest  that  inhibitory  synapses  might 
be  more  effective  in  their  activity  than  excitatory  synapses.  In  summary,  the  authors  establish  that 
‘a  random  walk  model  with  balanced  excitatory  and  inhibitory  inputs  allows  the  neuron  to  behave 
as  an  integrate  and  fire  device  and  maintain  a  reasonable  response  rate,  with  irregular  ISP. 

The  use  of  functional  magnetic  resonance  imaging  to  study  the  human  visual  cortex  directly  is  seen 
in  subsequent  papers.  In  [24],  Engel  et  al.  used  a  contrast  reversing  checkerboard  as  stimuli.  As 
the  stimulus  moved  from  fovea  to  periphery,  the  fMRI  signal  was  delayed  at  locations  containing 
neurons  with  peripheral  receptive  fields  when  compared  to  those  containing  neurons  with  foveal 
receptive  fields.  A  rotating  wedge  stimulus  was  used  to  study  the  retinotopic  organization  with 
respect  to  polar  angle.  They  also  obtained  clear  demarcation  between  various  areas  in  the  visual 
cortex.  Similar  stimuli  were  used  in  all  subsequent  fMRI  studies  by  other  authors  as  well.  A  study 
similar  to  [24]  of  the  human  visual  cortex  using  fMRI  was  conducted  by  Sereno  et  al.  in  [56].  The 
extent  and  borders  of  various  visual  areas  were  identified,  and  iso-eccentric  and  iso-polar  maps  of 
the  visual  areas  obtained.  A  comparison  between  human  and  other  primate  visual  areas  is  made, 
and  the  cortical  magnification  factors  plotted,  as  shown  in  Fig. 9.  Their  studies  imply  that  humans 
have  extreme  emphasis  on  the  center  of  gaze,  however,  they  leave  the  question  of  the  reason  for  the 
emphasis  (emphasis  on  the  retina  vs  emphasis  on  the  cortex)  open. 


Fig.  4.  Schematic  summa¬ 
ry  of  retinotopic  visual  ar¬ 
eas  in  the  owl  monkey  (7), 
the  macaque  monkey  [21), 
and  the  human  (present 
study)  at  the  same  scale 

(A)  and  approximately  nor¬ 
malized  by  the  area  of  VI 

(B) .  (Human  VI  is  twice  the 
area  of  macaque  VI ,  with 
larger  ocular  dominance 
columns  and  cytochrome 
oxidase  blobs,  but  a  similar 
number  of  cells.)  Visual  ar¬ 
eas  in  humans  show  a 
close  resemblance  to  visu¬ 
al  areas  originally  defined 
in  monkeys.  The  anterior 
border  of  the  visual  cortex 
in  humans  was  estimated 
by  using  the  superior  tem¬ 
poral  sulcus  and  intrapari- 
etal  sulcus  as  landmarks. 
In  (C),  the  mapping  func¬ 
tions  (heavy  lines;  scale  is 
on  the  left  axis)  and  magni¬ 
fication  factor  functions 
(light  lines;  scale  is  on  the 
right  axis)  are  shown  for 
the  upper  field  representa¬ 
tions  of  human  VI ,  V2,  VP, 
and  V4  (24).  The  VI  map¬ 
ping  functions  for  owl 
monkeys  (OM*,  dotted) 
and  macaque  monkeys 
(MM*,  dashed)  shown  at 
the  left  were  scaled  up  to 


Eccentricity  (degree) 


match  the  overall  size  of  human  VI  (25).  An  increased  emphasis  on  the  center-of-c 
evident.  CMF,  cortical  magnification  factor. 


}  in  human  VI  is 


Figure  B.9:  A  set  of  results  from  [56] 
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Balasubramanian  et  al.  [6]  modifies  the  previously  proposed  expressions  for  log  polar  mapping  with 
the  idea  of  a  dipole  mapping: 


w  =  log(z-fa)  -  log(z+b) 

where  w  is  a  point  on  the  cortical  plane,  z  is  a  point  on  the  visual  field  and  a  and  b  are  constants.  The 
authors  point  out  that  the  modified  expression  captures  the  shape  of  the  VI  boundary  exhibited  at 
peripheral  representation  as  well  as  the  fact  that  inverse  cortical  magnification  factor  is  sub-linear  in 
the  peripheral  field.  The  dipole  mapping  is  then  further  modified  in  order  to  explain  more  complex 
features  like  the  boundary  conditions  between  VI,  V2  and  V3  and  the  spacing  of  iso-eccentricity 
lines. 

Dougherty  et  al. [21]  describes  the  use  of  fMRI  to  study  visual  cortex  areas  VI,  V2  and  V3.  The 
authors  obtained  maps  of  the  locations  of  these  areas  in  the  cortex,  and  determined  their  surface 
areas  for  the  central  2°-12°  of  the  visual  field.  For  the  central  l°-2°,  the  various  areas  were  not 
distinguished,  and  the  combined  surface  area  for  this  central  region  of  the  visual  field  was  determined 
to  be  2100  mm2.  The  variation  of  cortical  magnification  with  eccentricity  was  found  to  be  in 
agreement  with  previous  literature.  The  plots  showing  the  variation  for  all  three  visual  areas,  for 
all  of  the  test  subjects,  were  obtained,  Fig.  10. 


Figure  4.  Surface  area  cortical  magnification  functions  for  all 
hemispheres,  (a)  VI .  (b)  V2  and  (c)  V3. 


Figure  B.10:  A  set  of  results  from  [21] 

Further  studies  of  the  human  visual  cortex  using  fMRI  are  seen  in  [50]  by  Schira  et  al.  The  authors 
determined  more  accurately  the  borders  between  VI,  V2  and  V3,  and  obtained  2  dimensional 
reconstruction  of  various  iso-polar  and  iso-eccentricity  lines  in  VI.  They  started  out  with  the 
following  expression  for  magnification: 


nr  =  k/  (r+a) 
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where  m  is  the  magnification  at  eccentricity  r,  and  k  and  a  are  constants.  They  determined  the 
value  of  a  to  be  approximately  0.75  for  both  VI  and  V2,  and  k  to  be  1.92  for  VI  and  1.59  for  V2. 
The  magnification  vs  eccentricity  plots  obtained  are  shown  in  Fig. 11.  The  authors  also  modified 
the  mapping  introduced  in  [6]  so  as  to  preserve  the  meridional  isotropy  of  the  area  at  the  expense 
of  small  degree  of  local  anisotropy.  The  question  of  whether  VI  is  actually  flat,  and  if  variability 
between  individuals  can  affect  the  construction  of  a  generic  model  for  VI  is  also  examined  in 
the  paper.  The  authors  make  the  case  for  VI  being  intrinsically  flat,  and  state  that  a  general 
model  is  the  best  way  to  estimate  VI  parameters,  in  the  cases  where  we  are  unable  to  measure  VI 
morphometry  for  individuals. 
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Figure  B.ll:  A  set  of  results  from  [50].  From  [50]  -  ‘Mean  vs. eccentricity  across  all  subjects 
fitted  by  Eq.l.  A:  measurements  across  areas  and  studies.  VI  data  are  plotted  in  red(fit:  a  = 
0.77  ±  0.03;k  =  19.2  ±  0.25  and  V2  data  are  plotted  in  blue(fit:  a  =  0.73  ±0.03;k  =  15.9  ±  0.21). 
The  gray  line  shows  the  cortical  magnification  function(Vl)  from  Horton  and  Hoyt  (1991a)(a  = 
0.75;k=17.6)  and  the  red  *  symbols  are  points  from  Dougherty  et  al.(2003),  error  bars  depict  SE. 
Because  our  measurements  did  not  extend  beyond  12°  of  eccentricity,  the  data  are  well  described 
by  the  monopole  model  or  the  corresponding  linear  magnification  function(Eq.l).B:  individual 
measurements  from  the  present  study.  Colored  graphs  show  the  results  in  individual  subjects,  the 
dotted  gray  line  the  fit  for  VI.  C:  mean  inverse  magnification  of  VI  with  the  best-fitting  straight 
line  for  ready  comparison  with  studies  such  as  Rovarno  and  Virsu  (1979,  Levi  et  al. (1984) ,  McKee 
and  Nakayama  (1984),  Stensaas  et  al.(2001),  Tyler  (2001),  and  Duncan  and  Boynton  (2003).’ 

Boucart  et  al. [15]  conducted  a  study  with  the  following  purposes  -  l)compare  the  performances  in 
implicit  and  explicit  recognition  tasks  using  same  set  of  stimuli  at  2  peripheral  eccentricities  -  30° 
and  50°  2)Assess  whether  people  who  have  lost  central  vision  early  in  life  develop  their  peripheral 
vision.  Implicit  recognition  involved  facilitation  for  old  stimuli  as  compared  to  new  stimuli  in  terms 
of  latency  and  accuracy.  Images  were  presented  in  2  stages  -  a  test  stage,  and  a  study  stage.  In 
the  study  stage,  the  images  presented  included  -  l)images  previously  presented  in  the  study  stage, 
2)images  of  objects  with  same  names  but  different  appearances  as  the  images  in  test  stage  -  for 
example,  the  image  of  a  different  cat  3)images  of  new  objects.  The  authors  found  no  significant 
difference  in  performance  between  identical  and  same  name  pictures  in  people  with  normal  vision. 
New  pictures  did  poorly  when  compared  to  identical  and  same  name  pictures  in  terms  of  response 
time  and  accuracy.  For  candidates  who  had  lost  central  vision,  identical  pictures  were  identified 
faster,  but  no  significant  difference  was  found  between  performances  for  new  pictures  and  same 
name  pictures.  In  people  with  central  vision  loss,  the  visual  system  appeared  to  have  put  its 
resources  into  developing  a  preferred  retinal  locus  which  would  stimulate  cortical  regions  which 
previously  responded  to  centrally  displayed  stimuli. 
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We  also  came  across  papers  on  contour  detection  during  the  course  of  our  investigation.  In  [27], 
Field  et  al.  present  their  results  from  5  experiments  which  tested  the  ability  of  observers  to  iden¬ 
tify  contours  from  within  a  field  of  randomly  oriented  elements.  The  experiments  were  conducted 
by  varying  orientation,  spacing  and  so  on  between  the  elements  belonging  to  the  contour.  Their 
results  led  the  authors  to  propose  a  local  association  field  in  our  visual  field.  In  [30],  Gintautas  et 
al.  measured  time  taken  by  human  observers  to  identify  segmented,  closed  contours  in  backgrounds 
with  clutter,  and  found  that  it  matched  their  model  of  local  association  fields  in  the  cortex.  When 
considering  probabilistic  methods  for  edge  detection,  it  has  been  found  that  there  is  a  local  maxi¬ 
mum  in  the  pairwise  probability  distribution  for  edge  elements  that  are  nearly  co-circular  as  well  as 
nearly  parallel.  The  experiments  in  [29]  suggest  that  humans  have  good  knowledge  of  the  pairwise 
statistics  of  edge  element  geometry  and  contrast  polarity,  and  that  they  use  the  knowledge  effi¬ 
ciently.  Feldman[25]  puts  forth  a  Bayesian  model  for  contour  detection,  and  suggests  that  humans 
detect  contours  using  such  a  model.  In  [45],  Ma  et  al.  suggest  that  humans  take  into  account  the 
reliability  of  the  observation  when  trying  to  detect  targets,  and  they  put  forward  a  log  likelihood 
ratio  to  capture  the  idea. 

B.4  Conclusion 

In  the  literature  we  came  across,  far  more  papers  support  the  idea  that  the  fovea  is  indeed  overrep¬ 
resented  in  the  cortex,  and  that  this  cannot  be  explained  by  the  higher  number  of  cones/ganglion 
cells  found  in  the  fovea  and  surrounding  areas  alone  [60] ,  [4] ,  [5] ,  [48] ,  [24] .  We  support  the  same 
hypothesis.  Fovea,  which  occupies  only  0.01%  -  0.025%  of  the  area  of  the  retina,  is  represented  by 
approximately  8%  of  the  striate  cortex  area,  and  this  emphasis  on  central  vision  does  not  directly 
correspond  with  the  emphasis  in  the  retinal  cell  density.  Hence,  the  additional  magnification  must 
happen  in  the  retino  cortical  pathway,  more  specifically  in  the  lateral  geniculate  nucleus  [5].  We 
put  forth  the  idea  that  the  additional  emphasis  for  central  vision  in  the  cortex  is  the  result  of  a 
learning  process  by  the  visual  system,  and  that  this  learning  was  triggered  by  the  higher  number  of 
cones/ganglion  cells  in  the  fovea  we  started  out  with.  Indirect  support  to  this  idea  can  be  seen  in 
[15].  In  the  absence  of  central  vision,  a  preferred  retinal  locus  which  triggered  the  original  central 
vision  cortical  regions  developed  in  the  test  subjects.  So  in  a  person  who  starts  out  with  normal 
central  vision,  the  fovea  and  surrounding  areas  would  get  such  a  preferential  treatment.  The  reason 
for  such  a  preference  in  the  first  place  would  be  the  higher  number  of  cones/ganglion  cells  in  the 
fovea,  which  naturally  makes  it  the  location  from  where  the  cortex  receives  maximum  information. 
A  preferential  treatment  would  lead  to  further  magnification  in  the  visual  pathway  and  the  cortex. 
We  have  completed  a  survey  of  the  important  pieces  of  literature  related  to  the  existence  of  a 
log  polar  transformation  in  the  human  visual  system.  The  idea  of  such  a  transformation  remains 
unchallenged,  and  researchers  continue  to  investigate  the  reasons  for  the  existence  of  such  a  trans¬ 
formation.  We  related  two  key  explanations  for  such  a  mapping,  and  looked  at  the  literature 
presenting  the  evidence  for  both  explanations.  We  support  one  of  the  two  explanations,  and  put 
forward  our  own  ideas  as  to  how  such  a  mapping  could  have  developed  in  the  visual  cortex. 
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Appendix  C 

Software 


C2LaTeX 

/*  This  is  version2.4  */ 


This  is  a  listing  of  the  source  code  for  the  simulator.  Included  are  some  support  functions  and 
debugging  functions  which  are  not  used  in  the  final  version. 


Created  by  Wesley  Snyder  on  3/20/12.  Copyright  (c)  2012  North  Carolina  State  University.  All 
rights  reserved.  This  program  simulates  area  VI  of  the  cortex  finding  a  straight  line  in  an  image. 

// 

SLDSimulation.c  Versionl.l  Created  June  18,  2012 
version2.0  Released  June  22,  2012 
Versionl.2  Created  June  23,  2012. 

Versionl.3  Created  June  25,  2012. 

Version  1.3:  This  version  reverses  the  order  of  searching  for  the  maxima  in  the  accumulator.  In 
this  version,  first,  the  accumulator  is  filtered  by  a  relatively  small  (3  x  3)  low  pass  filter.  Then,  an 
image  is  created  by  running  findpsr  find  peak  to  sidelobe  ratio.  Then  the  maxima  of  that  image  are 
located.  Versionl.4  Created  July  5,  2012. 

Version2.0  Created  July  9,  2012. 

Version  2.0  cleans  up  a  few  messy  bits  of  code,  speeds  up  a  few  things,  and  emphasizes  the  brightest 
peak  in  the  reconstructed  image. 

Version  2.1  allows  the  command  line  control  of: 

•  sigma  (of  the  edge  detector) , 

•  number  of  iterations 

•  some  other  stuff 

Version  2.2:  Created  July  18,  2012  Version  2.2  includes  most  of  the  functions  into  a  single  file 
(primarily  to  make  documentation  easy).  It  also  allows  increasing  the  Q  of  the  Gabor  filter  by 
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optionally  modifying  the  mu  loop  at  the  beginning  of  the  function  Accumulate.  The  option  is 
implemented  (at  this  release)  by  using  ifdef  Q 
Version  2,3,  Created  July  22,  2012 

This  version  finds  the  top  three  lines,  and  shows  them  in  the  plot. 
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C.l  Definitions  and  structures 


The  simulation  assumes  that  at  each  point  in  the  visual  field  there  are  p-^numorientations)  orienta¬ 
tions  covering  a  range  totaling  180  degrees.  This  is  in  agreement  with  literature  [51]  that  suggests 
the  orientation  sensitivity  of  higher  mammals  (including  presumably  humans)  is  about  ten  degrees. 
The  number  of  distinct  orientations  is  thus  18. 

The  program  goes  through  the  following  steps: 

1.  read  in  the  input  image  and  initialize  (function  initialize) 

2.  pass  the  input  image  through  a  temporal  high-pass  filter  (function  HighPass).  The  result  of 
this  is  that  over  time,  if  the  image  does  not  change  it  will  fade  away. 

3.  apply  a  spatial  high-pass  filter  to  the  image  (function  myflGabor). 

4.  For  every  pixel  in  the  resulting  high-pass  filtered  image,  calculate  its  contribution  to  an  ac¬ 
cumulator.  (function  Accumulate)  Note  that  there  are  seven  possible  modes  of  accumulation. 
Mode  5  seems  to  work  best,  and  is  the  default.  The  mode  is  set  by  using  the  -nr  switch  on 
the  command  line. 

5.  distinctive  points  in  the  accumulator  are  identified  and  marked  (function  findpsr).  This 
function  was  so  named  because  it  was  originally  designed  to  use  the  maxima  of  the  peak-to- 
sidelobe  ratio.  However,  it  turns  out  that  the  psr  doesn’t  work  very  well. 

6.  for  each  peak  in  the  accumulator  which  has  been  marked,  draw  a  line  on  the  screen. (function 
Houghinvsub) . 

The  accumulators  (there  are  two)  are  floating  point,  and  have  a  width  of  180  +  2P,  where  P  is  a  pad 
value  set  by  the  variable  PAD,  which  is  defined  int  Neuralhough.h.  The  accumulator  is  similarly 
padded  in  the  vertical  direction.  Currently  PAD  is  10  pixels.  The  purpose  of  the  padding  is  to  allow 
a  peak  on  occur  on  the  exact  left,  right,  top,  or  bottom  of  the  accumulator.  BY  CONVENTION, 
when  a  function  is  called  with  accumulator  coordinates  as  arguments,  the  actual  coordinates  will 
be  used  (rho  takes  on  values  between  —  p  and  p  and  theta  will  be  an  integer  between  0  and  180 


#def ine  MINBRIGHTNESS  50 
#def ine  LASTITERATION  30 


#include 

#include 

#include 

#include 


<f lip . h> 
<stdio . h> 
<math . h> 
<sys/stat .h> 


#include  <NeuralHough.h> 
struct  clp  { 

char  *infilename ; 
int  RecordFlag; 
int  mode ; 
int  loops; 
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float  stddev; 
int  degreespersample ; 
float  wavelength; 
float  aspect; 
int  numpeaks ; 


The  structure  param  is  used  to  hold  more  or  less  global  information  about  the  state  of  the  program, 
such  as  pointers  to  the  images,  before,  during,  and  after  high  pass  temporal  filtering,  edge  detector 
outputs,  etc. 
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C.2  Help  Function 


The  function  usage  explains  the  argument  list  to  the  program.  The  function  is  declared  static  to 
avoid  linking  name  conflicts. 

/*===================================================*/ 

/*  usage  */ 

/*===================================================*/ 

static  void  usage (void) 

{ 

printf ("Usage :  SLDSimulation  -i  filename  [— s]  [-v  var] [-m  x]  [-1  iterations]  -d  pixelspers; 
printf("  The  file  name  will  be  saved  as  retina. if s  if  retina. if s  changes,  \n"); 

printf ("  the  HT  will  be  updated. \n") ; 

printf ("  If  the  -s  flag  is  present,  it  will  indicate  saving  the\n"); 

printf ("  evoluation  of  the  accumulator  in  a  (rather  large)  file.\n"); 

printf ("  The  -m  switch  selects  the  mode  of  accumulation.  Default  is  5\n"); 

printf ("  The  -v  switch  allows  the  user  to  specify  the  std  dev  (sigma)  of  the  Gaussian 

printf ("  Default  is  0.75\n"); 

printf ("  The  -1  (loops)  switch  allows  the  user  to  specify  the  number  of  iterations.  F; 

printf ("  signal-to-noise  ratio,  so  the  system  is  most  sensitive  at  iteration  3. 

printf ("  the  -n  switch  controls  the  number  of  maxima  to  plot.  Default  is  3\n"); 

printf ("  The  -d  switch  controls  the  number  of  degrees/sample.  For  example,  since  the  ] 

printf ("  180  always,  degrees/sample  indirectly  also  selects  the  number  of  orien' 

printf ("  one  might  choose  degrees/sample  to  be  10,  which  would  cause  18  samples 

printf ("  per  samples  would  require  ten  samples.  Default  is  ten  degrees/sample .  1 

printf("  The  -w  switch  controls  the  wavelength  of  the  Gabor  filter.  Default  is  2.0.\n 

printf("  The  -a  switch  controls  the  aspect  ratio  of  the  Gabor  filter.  Default  is  0.75 

exit (-1) ; 

> 


65 


C.3  Parsing  the  Command  Line 


The  function  clparse  reads  the  command  line  and  interprets  the  command  arguments, 
initializes  several  variables 


/* - CLPARSE - */ 

/*  command  line  parser  */ 

static  void  clparse(int  argc,char  **argv, struct  clp  *stuff) 

{ 

int  i ; 

void  usage (void); 

/*  set  default  values  in  returned  structure  */ 
if(argc  ==  1) 

{ 

printf("You  must  specify  arguments  to  this  program\n"); 
printf("use  SLDSimulation  -h  for  more  details\n") ; 
exit (-1) ; 

> 

stuf f->inf ilename  = 

stuff ->RecordFlag  =0; 

stuff ->mode  =5; 

stuf f ->stddev=0 . 75 ; 

stuf f ->loops=3 ; 

stuf f->degreespersample  =  10; 

stuf f->wavelength  =  2.0; 

stuf f->aspect  =  .75; 

stuf f->numpeaks  =  3; 

ford  =  1;  i<  argc;i++) 

{ 

if  (* (argv [i] )  ! = 

{ 

printf ("Command  line  arguments  must  start  with  a  dash\n"); 
exit (-1) ; 

> 

switch  (  *(argv[i]  +1)) 

{ 

case  ’h’ : 

{ 

usage () ; 
exit (0) ; 

> 

break; 


It  also 
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case  1 i ’ : 

stuf f->inf ilename  =  argv[i+l]; 
i++; 

break;  /*  infilename  */ 
case  ’m’ : 

stuff->mode  =  atoi (argv  [i+1] ) ; 
i++; 

break;  /*  mode  */ 
case  ’s’  : 

stuf f->RecordFlag  =  1; 
break; 

case  ’1’: 

stuff->loops  =  atoi (argv [i+1] ) ; 
i++; 

break;  /*  loops  */ 
case  ’n’ : 

stuf f->numpeaks  =  atoi (argv  [i+1] ) ; 

if (stuf f->numpeaks  >6) 

{printf ("Cannot  draw  more  than  6  peaks. \n") 
exit (-1) ;> 

i++; 

break;  /*  number  of  peaks  to  be  shown  */ 
case  ’  v  ’  : 

stuf f->stddev  =  atof (argv [i+1] ) ; 

i++ ; 

break;  /*  stddev  of  edge  detector  blur  */ 
case  ’d’ : 

stuf f->degreespersample  =  atoi (argv  [i+1] ) ; 

i++ ; 

break;  /*  stddev  of  edge  detector  blur  */ 
case  ’w’ : 

stuf f->wavelength  =  atof (argv  [i+1] ) ; 

i++ ; 

break;  /*  Gabor  wavelength  */ 
case  ’a’ : 

stuf f->aspect  =  atof (argv [i+1] ) ; 

i++ ; 

break;  /*  Gabor  aspect  ratio  */ 
default : 

printf ("SLDSimulation:  unknown  option  \n") ; 

exit (-3) ; 

break; 

} 

> 

> 
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C.4  Saving  Accumulator  to  Disk 

This  is  a  function  used  primarily  for  debugging.  It  saves  the  floating  point  accumulator  to  disk  as 
a  floating  point  ifs  image.  The  padding  is  also  saved. 


/*=========================Sa.veAccPadded2Disk============================*/ 

/*  */ 
/*=================================================================*/ 

void  SaveAccPadded2Disk (float  **acc,char  ^filename , struct  param  *p) 

{ 

IFSIMG  hout; 
int  hr, he; 
float  value; 
int  hlen[3] ; 

void  copyandconvert (float  **, IFSIMG); 

hlen  [0] =2 ; hlen [1] =180+2*PAD ; hlen [2] =2*p->Accrows+2*PAD ; 
hout  =  ifscreate("f loat" ,hlen, IFS_CR_ALL,0) ; 

for(hr=0;hr  <  2*p->Accrows+2*PAD;hr++)f or (hc=0;hc  <  180+PAD;hc++) 

{ 

value  =  acc[hr]  [he]  ; 

//  if  (ace  [hr]  [he]  >  10 . 0)printf  ("writing  °/0f  at  °/0d  °/0d  to  padded  ace  file  \n" , value, hr  ,hc’ 

if sfpp(hout ,hr ,hc , value) ; 
if (isnan(if sfgp(hout ,hr ,hc) ) ) 

printf  ("SaveAccPadded2Disk:nan  at  °/0d  °/0d\n"  ,hr  ,hc) ;  ff lush  (stdout) ; 

} 

if spot (hout , filename) ; 
if sf ree (hout , IFS_FR_ALL) ; 
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C.5  Computing  a  histogram  of  the  Accumulator  values 


This  is  also  primarily  used  for  debugging.  It  prints  out  the  count  of  how  many  cells  in  the  accu¬ 
mulator  have  distinctive  brightness  as  a  function  of  angle. 


/*=================================HistogramAccuinulator===============*/ 

void  HistogramAccumulator (float  **Acc, struct  param  *p) 

{ 

int  r , c ; 
float  sum; 

for(c=0; c<180  +  2  *  PAD;c++) 

{ 

sum=0; 

for(r=0;r  <  2*p->Accrows  +  2  *  PAD;r++) 
sum  +=  Ac c  [r]  [c]  ; 

//  if  (sum  !=  744.0)  printf  (" iteration  %d,  °/0d  °/0f \n"  ,p->iteration,  c , sum)  ; 

} 

} 
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C.6  Main  Program 


The  program  uses  a  3D  data  set  to  store  the  output  of  the  orientation-sensitive  edge  detectors.  The 
name  of  that  image  is  Gaborout.ifs.  The  main  program  initializes  all  the  data  structures  and  runs 
the  program. 

The  input  image  is  read  from  a  disk  into  memory,  and  then  high-pass  filtered  with  respect  to  time. 
In  this  way,  if  the  same  stimulus  is  presented  only  once,  it  will  fade  over  time.  Each  time  through 
the  main  loop,  the  program  checks  the  creation  date  on  the  input  file  (retina. ifs).  If  the  file  has 
changed,  it  will  be  read  and  input  to  the  filter. 


/•  . ==== .  === .  ■  ■  .  === .  - .  *  / 

/*  */ 

/*  main  */ 

/*  */ 

/•  —  -  . —  -  .  - .  === .  .  »/ 


int  main (int  argc,  char  *argv[]) 

{ 


int  displayindex;  //  number  required  to  know  which  display  to  use 

struct  param  mparam,*p;  //  a  list  of  global  parameters 

p=&mparam; 

int  inititalize (char  *,  struct  param  *) ; 

int  HighPassFilter (struct  param  *,int); 

extern  void  copyandconvert (float  **,  IFSIMG) ; 

int  Gabors (IFSIMG  , IFSIMG,  float  , float, float) ; 

int  WriteToIFSDisplayWindow(int , IFSIMG, int , float , int , int) ; 

int  Accumulate (IFSIMG, float  **,  struct  param  * , int) ; 

void  writefile (IFSIMG  ,char  *  ,int  ); 

void  usage (void); 

int  findpsrl (struct  param  *) ; 

float  findpsr2 (struct  param  *) ; 

void  saveacc (struct  param  *) ; 
void  clparse(int  ,char  *[],  struct  clp  *) ; 


void  zap (IFSIMG); 

void  zapacc (float  **, struct  param  *) ; 
void  clparse (int , char  *[],  struct  clp  *) ; 
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void  clip (IFSIMG) ; 
struct  stat  mystat,  *info; 
int  oldmtime; 
int  status; 

struct  clp  myclp,  *stuff; 

The  accumulator  is  180  degrees  wide,  and  divided  into  NUMORIENTATIONS  blocks.  This  variable 
is  set  in  the  include  file  to  be  18.  However,  in  modifying  the  program,  be  careful.  The  numbers  180, 

18,  or  10  (degrees  per  sample)  may  be  used  without  reverence  to  this  macro.  This  inconsistency 
will  be  corrected  in  future  releases. 

It  is  important  to  remember  that  the  accumulator  is  actually  a  doubly-indexed  floating  point  array. 

This  array  has  dimensions  (2 PAD  +  180)  x  2 Accrows.  The  accumulator  itself  is  180  columns  wide, 
but  is  imbedded  in  a  larger  array  to  avoid  boundary  problems 

stuff =&myclp; 

clparse (argc , argv, stuff) ; 

//  printf ("SLDSimulation  version  2.1\n"); 

p->degreespersample  =  stuf f->degreespersample ; 
p->numpeaks  =  stuf f->numpeaks ; 

p->numorientations  =  180/stuf f->degreespersample ; 
inf o=&my stat ; 

p->df  =  180/stuf f->degreespersample ;  //  pixels  per  sample,  also  degrees  per  sample 
p->deltatheta  =  p->degreespersample  /  RAD2DEG;  //  space  betweeen  samples,  measured  in  radii 

oldmtime=0 ; 

The  first  argument  is  the  name  of  the  input  file.  That  is  stuff-/, infilename .  The  name  is  passed  to 
the  initialize  function.  First,  we  check  that  there  are  at  least  two  arguments. 

if (stuf f->RecordFlag)  p->RecordFlag  =  l;else  p->RecordFlag=0 ; 
p->iteration  =0; 

p->displayl  =  initialize (stuf f->inf ilename ,p) ; //  initalize  all  images 


The  input  will  be  read  if  it  needs  to  be  read  and  update  y.  Note  that  this  is  doing  high  pass  filtering 
in  time. 

//  printf ("Waiting  for  something  to  happen :\n"); 
do 
{ 

p->iteration++ ; 

//  printf  ("  °/0d\n"  ,p->iteration) ;  f  f  lush(stdout) ; 

stat (p->inf ilename , inf o) ; 
if  (inf o->st_mtime  >  oldmtime) 

{ 


/*  enter  here  if  the  input  file  has  been  rewritten  */ 
/*  read  the  new  input  image  */ 
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//  printf ("Something  Happened  at  iteration  °/0d!  ! "  ,p->iteration) ; 

oldmtime=inf o->st_mtime ; 
if sReRead(p->inf ilename ,p->retina) ; 
p->iteration=l ; 

} 

loopl:  status  =  HighPassFilter (p,p->iteration) ; 

if (status  ==0) 

{ 

p->iteration++ ; 
printf (" . ; 
goto  loopl; 

} 

function  Gabors  will  compute  the  Gabor  filter  of  the  input  at  all  the  orientations 


Gabors (p->y,p->v, stuff ->stddev, stuff ->wavelength, stuff ->aspect) ; 

zapacc (p->Accl ,p) ;  //  set  the  last  version  of  the  accumulator  back  to  zero 

Here  is  the  call  to  the  critical  function:  Accumulate.  The  function  accumulate  will  add  up  all  the 
inputs  to  the  accumlator 


Accumulate (p->v,p->Accl ,p, stuff ->mode) ; 

the  following  two  lines  have  been  commented  out  to  make  the  program  run  faster,  but  are  very 
useful  for  debuggin. 


//  HistogramAccumulator (p->Accl ,p) ; 

//  SaveAccPadded2Disk(p->Accl , "HoughAfterAccumulate . if s" ,p) ; 


the  function  saveacc  is  only  called  for  debugging.  It  writes  a  large  file.  Before  opening  the  output 
file,  the  function  saveacc  will  check  to  see  if  the  user  wants  to  write  the  file.  This  desire  is  expressed 
by  setting  the  -s  switch 


saveacc (p) ; 

copyandconvert (p->Accl ,p->iAcc) ; 

//  if spot (p->iAcc , "iAcc. ifs") ; 

here  is  where  the  accumulator  is  displayed  on  the  screen,  at  the  upper  left  of  the  screen 
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Wr it eTo IFSDi splayWindow (p->di splay 1 ,p->iAcc ,0 , 1 . 0,0 , 0) ; 
sleep(l) ; 

//  usleep(250000) ;  //  wait  250  microseconds 


now  find  the  point  with  maximum  peak-to-sidelobe  ratio  the  function  findpsr  finds  the  peaks  of  an 
accumulator  at  a  particular  point 


if (  findpsrl(p)  >  0)  showHoughinv(p) ; 

}  while  (p->iteration  <  stuf f->loops) ; //  end  while  loop 
sleep(5) ; 

}//  end  main 
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C.7  saveacc 


/*===================saveacc==========================*/ 

/*  writes  a  series  of  accumlators  to  a  3D  float  disk  file  */ 

/*  */ 
void  saveacc (struct  param  *p) 

{ 

double  v; 
int  r , c ; 

float  max, min, scale ; 

int  len[3] ; 

char  filename [132] ; 

IFSIMG  outimg; 

if (p->RecordFlag  ==  0  ){/*printf ("Returnflag  =  0") ; */return; } 
#if def  PNGS 

len  [0] =2 ; len  [1] =180 ; len [2] =2*p->Accrows ; 
outimg=ifscreate("u8bit" , len, IFS_CR_ALL , 0) ; 
for(r=0;r  <  2*p->Accrows ; r++) 
for(c=0;c  <  180;c++) 

{ 

v=p->Accl [r+PAD]  [c+PAD] ; 
if sfpp( outimg, r ,c,v) ; 

} 

sprintf  (filename ,  "saveacc°/0d.png"  ,p->iteration) ; 
if sCVpot( outimg, filename) ; 
if sfree (outimg, IFS_FR_ALL) ; 

#else 

//  printf ("Entering  saveacc\n" ) ; f f lush(stdout) ; 
if (p->iteration  ==  LASTITERATIQN) 

{ 

//  printf ("saveacc  saving\n") ; f f lush(stdout) ; 

//  if spot (p->savedacc , "savedacc . ifs") ; 

return; 

} 

//  printf ("saveacc  adding  a  frame\n") ; f f lush(stdout) ; 

//  find  max  and  min  in  this  frame 
max  =  -100000 . 0 ;min=100000 . 0 ; 
for(r=PAD;r  <  2*p->Accrows+PAD ; r++) 
for(c=PAD;c  <  180+PAD;c++) 

{ 

v  =  p->Accl[r]  [c] ; 
if (v>max)  max  =v; 
if (v<min)  min  =v; 

> 

if (max  !=  min)  scale  =  255.0  /  (max  -  min); else  scale  =1.0; 
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//  now  scale  the  image  as  it  is  stored. 


for(r=PAD;r  <  2*p->Accrows+PAD ; r++) 
for(c=PAD;c  <  180+PAD;c++) 

{ 

v=p->Accl[r]  [c] ; 

ifsfpp3d(p->savedacc,p->iteration,r-PAD,c-PAD, (v-min) *scale) ; 

//  if (r  ==  343  &&  c  ==  125) 

//  printf  ("saveacc:°/0d  °/„d  °/0d  °/of\n"  ,p->iteration,r,c,v) ; f f lush(stdout) ; 

> 

#endif 

//  printf ("Leaving  saveacc\n") ; f f lush(stdout) ; 

> 
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C.8  Find  peaks  in  accumulator 


/*======================f indpsr 1=======================  */ 

this  function  searches  the  accumulator  to  find  maxima  and  if  */  that  is  over  a  threhold,  marks  it 
in  the  accumulator  and  erases  all  other  points 
marking  is  done  by  making  the  point  negative 

/*=======================================================*/ 

#def ine  PSRTHRESH  1.0  //  no  longer  used 
int  f indpsr 1 (struct  param  *p) 

{ 

float  maxpsr ,pssr ; 
int  maxpsrrow,maxpsrcol ; 
int  r,c, count; 
float  value; 

extern  float  psr (float  ** , int , int , int , int) ; 
extern  float  psraccl (struct  param  *p,int  r,int  c  ); 
extern  void  remarkaccl (struct  param  *) ; 

extern  void  ShowAccNeighborhood (struct  param  *, float  **,int,  int); 
extern  int  localmax (struct  param  *, float  **,int  ,int); 
void  copyandconvert (float  **,IFSIMG); 
count  =  0; 

first  copy  all  of  Accl  to  Acc2 

for(r=0;  r  <  2*p->Accrows+2*PAD  ;r++)  for(c=0;  c  <  180+2+PAD ; C++) 
p->Acc2[r]  [c]  =p->Accl  [r]  [c]  ; 

Now,  Acc2  is  the  acc.  Search  Acc2  for  maxima.  When  a  max  is  found,  the  corresponding  in  Accl 
is  made  negative. 

for(r=-p->Accrows;  r  <  p->Accrows  ;r++)  for(c=0;  c  <  180;c++) 

{ 

if(  p->Acc2 [r+p->Accrows+PAD]  [c+PAD] >100 . 0  &&  localmax(p,p->Acc2,r,c)  ) 

{ 

//  if (c  ==  145) 

//  { 

//  printf ("Acc2 : ") ; 

//  ShowAccNeighborhood (p, p->Acc2, r , c) ; 

//  printf ("Accl :") ; 

//  ShowAccNeighborhood(p,p->Accl,r,c) ; 

//  > 

value  =  p->Acc2  [r+PAD  +  p->Accrows]  [c+PAD] ; 
if (value  >  10000.0  &&  localmax(p,p->Acc2,r,c)) 
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{ 

p->Accl [r+PAD  +  p->Accrows]  [c+PAD]  = 
value  *  -1.0; 

//  printf  ("f ind:  marking  accl=%f  at  rho  =  °/0d  theta  =  °/od,  (°/0d  °/0d)\n", 

//  value,  r, c,r+p->Accrows+PAD, c+PAD) ; f f lush(stdout) ; 

count++ ; 

} 

> 

else 

{ 

p->Accl [r+PAD  +  p->Accrows]  [c+PAD]  = 
p->Acc2  [r+PAD  +  p->Accrows]  [c+PAD] ; 

} 

> 

//  printf  ("findpsrl :  marked  °/0d  points\n" , count) ; fflush(stdout) ; 

//  remarkaccl (p) ; 

#def ine  SHOWINPUTACC 
#if def  SHOWINPUTACC 
{ 

void  SaveAccPadded2Disk(f loat  **acc,char  ^filename , struct  param  *p) ; 
SaveAccPadded2Disk(p->Acc2 , "HoughPadded . if s" ,p) ; 

> 

#endif 

return  count ; 

> 

/*======================f indpsr2 

/*  not  used  in  this  version 
/*============================== 

float  findpsr2 (struct  param  *p) 

{ 

float  maxpsr, pssr; 
int  maxpsrrow,maxpsrcol ; 
int  r , c ; 

extern  float  psr (float  ** , int , int , int , int) ; 
printf ("Entering  findpsr2\n") ; f f lush(stdout) ; 

maxpsr  =  -100000.0; 

for(r=PAD;  r  <  2*p->Accrows  +PAD  ;r++)  for(c=PAD;  c  <  180+PAD;c++) 

{ 

pssr=p->Accl [r] [c] ; 
if(pssr  >  maxpsr) 

{ 

maxpsr  =  pssr; 


=*/ 

*/ 

*/ 
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maxpsrrow  =  r; 
maxpsrcol  =  c; 

> 

> 

printf  ("Maximum  Acc  value  of  °/0f  found  at  row  =  °/0d  (rho  =  °/0d)  theta  =  °/0d\n" , 
maxpsr , 

maxpsrrow,  maxpsrrow- (p->Accrows) -PAD ,maxpsrcol-PAD) ;  f f lush(stdout) 
return  maxpsr; 

> 
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C.9  printmatrix 

This  is  a  debugging  utility,  used  only  in  this  program 


static  void  printmatrix (double  **m,int  nr,int  nc) 

{ 

int  row, col; 

for (row  =  l;row  <=  nr;row++) 

{ 

f or (col=l ; col  <=nc;col++) 

printf("°/0f  "  ,m[row]  [col] ) ; 
printf ("\n") ; 

} 

> 


79 


C.10  Initialization 


The  initialize  function  sets  things  up 


/*============================================================*/ 

/*  initialize  */ 

/*============================================================*/ 


The  initialize  function  sets  up  everything  that  is  required  it  reads  the  input  image  file  and  computes 
an  accumulator 

int  initialize (char  ^filename , struct  param  *p) 

{ 

IFSIMG  input; 
int  len[4] ; 
int  displayindex; 

int  CreateIFSDisplayWindow( IFSIMG, float, int, int) ; 
extern  void  copyandconvert (float  **,  IFSIMG); 
int  ifsany2any(IFSIMG, IFSIMG) ; 
int  f ; 

float  **  SetUpAccptr (IFSIMG) ; 
float  **MakeAccumulator (int) ; 

extern  void  testaccindexing(struct  param  *) ; 
float  makekernel (float  ^kernel , int  kernelradius) ; 

//  printf ("Entering  Initializa\n") ; f f lush(stdout) ; 


Compute  the  f3  matrix  which  is  used  to  estimate  the  position  of  the  peak  more  accurately  in  the 
accumulator,  this  is  accomplished  as  follows:  Assume  there  are  n  possible  orientations  which  can  be 
estimated  by  an  orientation-sensitive  Gabor  function.  Further,  assume  the  directions  of  maximum 
sensitivity  of  each  of  those  Gabor  functions  are  defined  by  direction  vectors  b\,  62, ...,  bn  is  projected 
onto  each  of  these  direction  vectors  to  produce  a  projection  fi, .  The  unknown  gradient  vector  G  is 
projected  onto  each  direction  vector  obeying 

m  =  bfG 

Collect  all  the  direction  vectors  into  a  single  matrix  B  =  [b±,  62,  •  ••,  bn\.  Then  the  projection  of  all 
the  points  can  be  collected  into  a  single  vector  equation 

BtG  =  n 

.  Which  can  be  solved  for  a  minimum-squared-error  estimate  of  G  by  (showing  matrix  sizes) 

£*2x1  =  ( BB  )2x2"®2xn/^nx  1 
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In  the  following  code,  the  matrix  ( BBT )  1B  =  f3  is  computed  and  will  be  stored  in  the  structure  p 
This  implementation  has  not  been  recently  tested.  It  should  work. 


{ 

double  **B , **Btrans , **BBtrans , **BBtransinv; 

dmatrix  is  rows,  columns,  in  that  order 

B=dmatrix(l ,2, 1 ,p->numorientations) ; 

Btrans=dmatrix(l ,p->numorientations ,1,2) ; 

BBtrans=  dmatrix (1 ,2, 1 ,2) ; 

BBtransinv=dmatrix(l , 2,1,2); 

p->beta  =  dmatrix (1 , 2 , 1 ,p->numorientations) ; 

There  may  be  some  confusion  with  indices  below.  I  use  matrices  indexed  l...n,  whereas  in  the  other 
parts  of  the  program,  indices  range  from  0  ...  n-1. 

for(f  =  0;  f  <  p->numorientations ; f ++) 

{ 

B[l][f+1]  =  cos (p->deltatheta*  f ) ; 

B[2][f+1]  =  sin(p->deltatheta*  f ) ; 


> 

#if def  TESTB1 

B  [1]  [1]  =1; 
B  [2]  [1]  =0; 
B [1] [2]  =0; 
B [2]  [2]  =1; 
B [1]  [3]  =0; 
B  [2]  [3]  =0; 
B [1]  [4]  =0; 
B  [2]  [4]  =0; 

B [1] [5]  =0; 
B  [2]  [5]  =0; 

B [1]  [6]  =0; 
B  [2]  [6]  =0; 

B [1]  [7]  =0; 
B  [2]  [7]  =0; 

B [1]  [8]  =0; 
B  [2]  [8]  =0; 

B [1]  [9]  =0; 
B  [2]  [9]  =0; 
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B  [1]  [10]  =0 
B [2]  [10]  =0 
B [1]  [11]  =0 
B  [2]  [11]  =0 

B  [1]  [12]  =0 
B  [2]  [12]  =0 
B  [1]  [13]  =0 
B  [2]  [13]  =0 
B  [1]  [14]  =0 
B [2]  [14]  =0 

B  [1]  [15]  =0 
B  [2]  [15]  =0 
B  [1]  [16]  =0 
B  [2]  [16]  =0 
B  [1]  [17]  =0 
B  [2]  [17]  =0 
B  [1]  [18]  =0 
B  [2]  [18]  =0 


#endif 

transpose (B , 2 ,p->numorientations , Btrans) ; 

#if def  TESTB 

printf("The  matrix  Btrans  is\n") ; f f lush(stdout) ; 
printmatrix (Btrans ,p->numorientations , 2) ; 

#endif 

if smatmult (B , Btrans ,BBtrans ,2 ,p->numorientations ,p->numorientations , 2) ; 
#if def  TESTB 

printf("THe  matrix  BBtrans\n") ; f f lush(stdout) ; 
printmatrix (BBtrans, 2, 2) ; 

#endif 

if sinverse (BBtrans ,BBtransinv, 2) ; 

#if def  TESTB 

printf ("BBransinv  is\n") ; f f lush(stdout) ; 
printmatrix (BBtransinv, 2 , 2) ; 

#endif 

if  smatmult (BBtransinv , B , p->beta , 2,2,2, p->numorient at ions ) ; 

#if def  TESTB 

printf ("beta  is\n") ; f f lush(stdout) ; 
printmatrix (p->beta,  2,p->numorientations) ; 

#endif 

} 
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p->retina=if spin(f ilename) ; 
p->inf ilename=f ilename ; 
if (p->retina  ==  0) 

{ 

printf  ("Could  not  read  input  file  °/0s\n" , filename) ; 
exit (-1) ; 

} 

if (p->retina->if sdims  !=  2) 

{ 

printf ("\nThis  program  can  only  process  two-dimensional  images\n"); 
printf ("You  can  extract  a  2D  frame  from  a  3D  image  using  ChooseFrame\n") ; 
exit (0) ; 

} 

len[0]  =2; 

len[2]=p->nr  =  if sdimen(p->retina, 1) ; 
len[l]=p->nc  =  if sdimen(p->retina, 0) ; 
if (p->retina->dtype  !=  IFST_32FLT) 

{ 

int  status; 

p->xnml  =  if screate("float" , len, IFS_CR_ALL , 0) ; 

status=if sany2any (p->retina,p->xnml) ;  //  create  a  floating  version  of  the  input 
//  printf  ("any 2 any  returned  °/0d\n" ,  status)  ; f f  lush(stdout) ; 

} 

else 

{ 

p->xnml  =  if screate("float" , len, IFS_CR_ALL , 0) ; 

} 

//  printf  ("Init :  0 . 2 ,  len=°/„d  °/„d  °/0d\n" ,  len  [0]  ,  len  [1]  ,  len  [2] ) ;  f  f  lush(stdout) ; 

The  output  file  is  set  up  in  such  a  way  that  the  data  will  not  be  compressed.  Thus,  if  it  must  be 
read  multiple  times,  it  will  be  fater  to  ,  it  will  be  fast 


p->retina->comp=0 ; 

if spot (p->retina, "retina. if s" ) ;  //  save  the  uncompressed  retina  image 

//  create  the  rest  of  the  images 

p->xn=if screate("float" ,len,IFS_CR_ALL,0) ; 
p->y=if screate("float" ,len, IFS_CR_ALL,0) ; 
p->iy=if screate("u8bit" ,len,IFS_CR_ALL,0) ; 
p->templ=ifscreate("float" ,len,IFS_CR_ALL,0) ; 
p->temp2=ifscreate("float" ,len,IFS_CR_ALL,0) ; 
len[0]  =3; 

len  [3] =p->numorientations ; 

//  create  a  convolution  kernel  to  be  applied  later  to  blur  the  accumulator 
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int  eye, jay; 
int  inducks ; 

p->kernelradius  =  1;//  a  radius  of  1  produces  a  3x3 

p->kernel  =  (float  *)malloc ( (p->kernelradius  *2+1)  *  (p->kernelradius  *  2  +1)  *  sizeof 
makekernel (p->kernel , p->kernelradius) ; 

printf("the  radius  °/0d  convolution  kernel  is  :  \n"  ,p->kernelradius) ; 
inducks=0 ; 

f or (eye=0 ; eye<l+2*p->kernelradius ; eye++) 

{ 

for(jay=0; jay<l+2*p->kernelradius; jay++)  printf ("%f  " ,p->kernel  [inducks++] ) ; 
printf ("\n") ; 

> 

} 

f f lush(stdout) ; 

p->v=if screate("float" ,len, IFS_CR_ALL,0) ; 


create  the  accumulator 


first,  figure  out  how  many  rho  values  there  are  Accrows  does  NOT  include  padding 


p->Accrows=  sqrt(len[l]  *  len[l]  +  len  [2]*  len[2]); 

//  printf  ("Initializing  Accrows  to  °/0d\n"  ,p->Accrows)  ; 
len  [0] =2 ; 
len  [1] =180 ; 

len [2] =p->Accrows  *  2;//  allocation  for  negative  rho  and  padding  will  be  done  in  MakeAccumi 
p->Acccols  =  len[l] ; 

p->Accl=MakeAccumulator (p->Accrows) ; //  accrows  does  not  include  the  doubling  of  vertical  s: 
p->Acc2=MakeAccumulator (p->Accrows) ; //  that  is  added  in  the  make 


Construct  an  ifs  image  which  is  big  enough  to  hold  the  entire  accumulator  including  the  padding 

len [0] =2 ; len [1] =180  +  2  *  PAD;  len [2]  =  2  *  p->Accrows  +  2  *  PAD; 
p->iAcc=ifscreate("u8bit" , len, IFS_CR_ALL,0) ; 


We  will  want  to  display  the  image,  reconstructed  from  the  Hough  so  we  will  create  an  image  to 
store  those  results  First,  figure  out  how  many  rows.  Thats  easy  we  will  use  the  number  of  rows  in 
the  accumulator.  Thats  more  than  necessary,  but  who  cares? 


len  [2] =p->Accrows  *  M_SQRT2  +  1.0; 
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Lets  make  the  number  of  columns  the  same  number 


len[l]  =  len[2]  ; 

p->reconstructed=if  screate ( "u8bit " , len , IFS_CR_ALL , 0) ; 
p->tempacc=if screate("f loat" ,len, IFS_CR_ALL,0) ; 

copyandconvert (p->Accl ,p->iAcc) ; 

p->displayl  =  CreateIFSDisplayWindow(p->iAcc , 1 . 0 , 0 , 0) ;  //  set  up  display  of  Acc 

The  second  display  is  used  to  display  the  reconstructed  lines 

p->display2  =  CreateIFSDisplayWindow(p->reconstructed, 1.0, 300,0) ; 

the  RecordFlag  used  below  identifies  whether  or  not  to  save  the  accumulator  to  a  (rather  large) 
diskfile 

if (p->RecordFlag  ==  1) 

{ 

len [0] =3 ; len [1]  =180 ; len [2] =2*p->Accrows ; len [3] =LASTITERATI0N-1 ; 
p->savedacc  =  if screate("float" ,len, IFS_CR_ALL,0) ; 

} 

//  p->Yimg  =  CreateIFSDisplayWindow(p->retina, 1 . 0) ;  //  set  up  display  of  Acc 
//  printf ("Returned  from  create\n") ; f f lush(stdout) ; 

//  test  accumulator 
//  testaccindexing(p) ; 

return  p->displayl; 

}  //  done  with  function  initialize 
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C.ll  Set  up  connections  of  Axons 


An  array  of  axons  is  modeled  by  a  triply-indexed  array  of  pointers  to  accumulator  cells.  Each 
pointer  models  a  single  axon.  The  three  indices  are  row,  column,  and  orientation.  This  function  is 
not  implemented  at  this  release 

#if def  AXONARRAY 

/*===================SetUpAxonArray========================*/ 

/*  create  an  array  of  pointers  to  the  accumulator  */ 

SetUpAxonArray (struct  param  *p) 

{ 

float  theta, dtheta, cth, sth; 
int  f ; 

dtheta  =  M_PI  /  (p->numorientations  *  p->df ) ;  //  how  many  frames  is  also  how  many  angles 
for(f=0;f  <  p->numorientations ; f ++) 

{ 

theta=f *dtheta; 

cth  =  cos(theta);  sth  =  sin(theta) ; 
f or (r=p->nr ; r>=0 ; r — )f or (c=p->nc ; c>=0 ; c — ) 

{ 

pointer  =  p->v->ifsptr  + 

> 

} 

> 

#endif 


86 


Figure  C.l:  The  input  image  is  high-pass  filtered  producing  an  output  image  y,  which  fades,  the 
decay  rate  depends  on  the  sampling  time  and  the  constant  b\.  Here,  a0  =  1,  a\  =  —1,  and  b\  =  .95 

C.12  High  Pass  Filtering 

C.12.1  Dynamic  Response  of  the  input  system 

The  dynamic  response  of  the  visual  system  is  nicely  summarized  by  Martinez-Conde  et  al.  [47]: 
“Even  our  own  visual  system  can  detect  stationary  objects  only  because  the  image  projected  onto 
our  retinas  are  never  stationary  for  long.”  That  a  stationary  image  fades  is  a  phenomenon  referred 
to  as  “visual  fading.”  [49,  20].  The  fact  that  scenes  we  observe  do  not  fade  from  view  is  commonly 
attributed  to  either  head  motion  of  microsaccadic  motion  of  the  eyes  themselves. 

In  this  project,  we  do  not  attempt  to  model  microsaccadic  motion,  and  therefore  must  allow  constant 
images  to  fade.  In  the  current  version  of  the  program,  we  also  are  not  using  camera  input,  but 
instead  use  a  single  file  as  input.  The  creation  time  of  that  file,  named  “img.ifs,”  is  polled  each 
iteration  of  the  program,  and  the  file  is  read  only  if  it  has  changed.  If  change  has  occurred,  it  is 
read  into  an  input  buffer  named  x. 

Whether  read  from  disk  or  not,  the  input  image  x  is  subjected  to  a  simple  recursive  high  pass  filter 

yn  =  a0xn  -  aix71”1  +  blVn~l  ,  (C.l) 

producing  output  y,  as  illustrated  in  Figure  C.l.  The  input  processing,  and  the  formation  of  the 
fading  by  high  pass  filter  is  performed  in  blocks  I  and  II  of  the  flow  chart  in  Figure  ??. 

In  block  I  of  that  figure,  the  assignment  xn-\  <—  xn  is  denoted  xnml  :==  xn,  and  is  implemented 

using  almost  no  computing  resources  by  simply  observing  that  these  are  not  images  but  pointers 

to  images,  and  swapping  them  in  the  following  way: 

temp  =  xnml; 

xnml  =  xn; 

xn  =  temp; 


/* - == - - - */ 

/*  HighPassFilter  */ 
/*  returns  0  if  the  output  image  is  null  */ 
/* - -= - */ 


int  HighPassFilter (struct  param  *p,int  iteration) 

{ 


void  zap(IFSIMG); 

int  if nullimage(IFSIMG, float) ; 

switch (iteration) 
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{ 

case  1 : 

zap(p->y) ;zap(p->xnml) ;flcp(p->retina,p->xn) ;  //  initialization 
return  0; 
break; 
case  2: 

f lcp(p->xn,p->y) ; 
f lcp(p->xn,p->xnml) ; 
f lcp(p->retina,p->xn) ; 
if (ifnullimage(p->y , 0) )  return  0; 
break; 
default : 

f lmults (p->y ,p->templ , 0 . 9) ;  //  fade  by  this  much 
f lcp(p->templ ,p->y) ; 

if (ifnullimage(p->y , 0) )  {printf ("def alt \n") ; return  0;} 
break; 

} 

return  1  ; 

> 
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C.13  Non  Maximum  Suppression 

function  to  suppress  edge  responses  which  are  not  maxima 


/  •  - . - .  -  . --  . •  / 

/*  nms  */ 

/• .  ■  ■  ■  .  ■  -  ■  -  .  . •  / 


void  nms(IFSIMG  inimg, float  angle) 

{ 

int  r , c ; 

float  rl,cl,r2,c2; 

int  nr,nc,irl,icl,ir2,ic2; 

float  d,v,dmax,max; 

extern  double  mysin(int) ,  mycos(int); 

nr  =  if sdimen(inimg, 1) ;nc  =  if sdimen( inimg, 0) ; 
for(r  =  3;r<  nr-3;r++)f or (c  =  3;c<  nc-3;c++) 

{ 

v  =  if sfgp( inimg, r,c) ; 
if (v  >  20.0) 

{ 

max  =  0.0; 

for(d=-2.0;d  <=  3.0;d+=1.0) 

{ 

rl  =  r+  d*  sin(angle) ; cl=c  +  d*cos (angle) ; 
irl=rl ; icl=cl ; 
v=if sfgp(inimg, irl , icl) ; 
if (v  >  max) 

{ 

max=v ; 
dmax  =  d; 

} 

} 

for(d=-2.0;d  <=  3.0;d+=1.0)  if (d  !=  dmax) 

{ 

r2  =  r+  d*  sin(angle) ; c2=c  +  d*cos (angle) ; 
ir2=r2 ; ic2=c2 ; 

//  printf  ("nms :  suppressing  °/0d  °/„d  because  °/0f  <  %f  \n" ,  ir2 ,  ic2 ,  if  sfgp(inimg,r2 ,  c2) 

ifsfpp(inimg,r2,c2,0.0) ; 
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> 


} 
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C.14  Calling  the  Gabor  Edge  Detectors 


The  Gabor  filter,  applied  to  a  point  xy  in  an  image,  computes 


g(x,y;  A,  (9,  ^,<7, 7)  =  exp  - 


x'2  +  7y2 


2a2 


cos  (  27 r—  +  i/j 
A 


Here,  the  Gabor  is  implemented  by  using  the  ifs  function  flGabor  The  normal  to  the  line  has  the 
angle  theta,  not  the  line  itself  A  2D  Gabor  kernel  implemented  by  flGabor  is  mathematically 


defined  as:  G(x,  y )  =  exp 
variable  aspect  where 


cos  (27rA)  where  A  is  the  variable  wavelength  and  7  is  the 


(C.2) 

(C.3) 


x'  =  x  cos  6  +  y  sin  6 
y  =  — x  sin  6  +  y  cos  6 


The  parameters  involved  in  the  construction  of  a  2D  Gabor  filter  are: 

•  The  variance  of  the  gaussian  function 

•  The  wavelength  of  the  sinusoidal  function 

•  The  orientation  of  the  normal  to  the  parallel  stripes  of  the  Gabor  function 

•  The  spatial  aspect  ratio  specifies  the  ellipticity  of  the  support  of  the  Gabor  function.  For  , 
the  support  is  circular.  For  the  support  is  elongated  in  the  orientation  of  the  parallel  stripes 
of  the  function. 

Note  that  flGabor  does  not  give  the  user  the  options  of  selecting  t/>,  and  therefore  this  version  of 
flGabor  is  centered  around  the  center  (since  it  is  a  cosine  function.  It  will  therefore  select  lines,  not 
edges. 


/  •  -  -  . - . - .  . - . -  -  . --•  / 

/*  Gabors  */ 

/* . -  - .  - .  .  -  -  -  -•/ 


int  Gabors(IFSIMG  y,IFSIMG  v,  float  stddev, float  wavelength, float  aspect) 

{ 

int  i;  //counter 
int  len[3] ; 

int  nf,nr,nc;  //number  of  output  frames 

int  nfloats;  //number  of  floats  in  a  single  frame 

float  *oldifsptr;  //save  the  old  ifsptr 

float  *fptr;  //  temportary  pointer  to  a  float 

double  dtheta;  //  increment  theta 

int  frame; 
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extern  int  myflGabor(IFSIMG,IFSIMG, float, float, float, float) ; 
char  newname  [32] ; 
len[0]  =2; 

nf  =  if sdimen(v,2) ;  //  get  number  of  frames  in 
len [1] =nc=if sdimen(v, 0) ;  len [2] =nr=if sdimen(v, 1) ; 
nfloats=  nr*nc; 
dtheta  =  M_PI/nf; 

oldif sptr=(f loat  *)v->ifsptr;  //remember  the  old  ifsptr 
fptr  =oldifsptr; 
f or (i=0 ; i<nf ; i++) 

{ 

float  angle; 
angle  =  dtheta*i; 
v->if sdims=2 ; 

v->ifsptr  =  (char  *)fptr; 

The  first  argument  is  the  input  ifs  image.  The  second  argument  is  the  3d  Gabor  output  image.  The 
output  image  has  180/stuff-^degreespersample  frames,  on  for  each  angle.  The  third  argument  is  the 
standard  deviation  of  the  Gaussian  of  the  filter,  the  fourth  argument  is  an  angle  of  the  gradient, 
in  radians.  An  angle  of  0  will  be  maximally  sensitive  to  vertical  edges.  The  final  argument  is  not 
used,  but  is  retained  for  compatibility  with  the  ifs  flGabor  function. 

#if def  EDGEGABOR 

myflGabor(y,v,stddev, dtheta  *  i,  wavelength); 

#else 

f lGabor(y,v, stddev, dtheta  *  i,  wavelength, aspect) ; 

#endif 

//  f labsolute (v, v) ;  //  take  magnitude  of  Gabor 

//  nms (v, dtheta*i) ;  //  do  nonmaximum  suppression  on  the  edge  image 

fptr  +=  nfloats  ; 

} 

//restore  the  3d  nature  of  the  image 
v->ifsptr  =  (unsigned  char  *) oldif sptr; 
v->if sdims=3 ; 
ifspot(v, "Gaborout .ifs") ; 

//  debug  it 
#if def  DEBUG1 

for (frame=0; frame  <  p->numorientations ; f rame++) 

{ 

float  max; int  rmax,cmax,r,c; 

printf  ("looking  at  frame  °/0d\n" , frame) ; 

max  =  -10000.0; 

for(r=0;r  <  nr ;r++)f or (c=0 ; c<nc ; C++) 
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{ 

//  printf ("Reading  °/„d  %d\n",r,c); 

if (if sf gp3d(v,f rame ,r , c)  >  max) 

{ 

max  =  ifsfgp3d(v,frame,r,c) ; 
rmax  =  r ; cmax  =c ; 

} 

} 

printf  ("Max  at  row  °/0d  col  °/0d  in  frame  °/od\n" , rmax,  cmax,  frame ) ; 
printf  ("frame  °/„d  row  %d  col  °/0d  =  °/0f  °/0f  °/0f  %f  °/0f\n"  .frame , rmax,  cmax, 
if sfgp3d(v, frame, rmax, cmax  -2), 
if sfgp3d(v, frame, rmax, cmax  -1), 
if sfgp3d(v, frame, rmax, cmax) , 
if sfgp3d(v, frame, rmax, cmax  +1), 
ifsfgp3d(v, frame, rmax, cmax  +2) 

) ; f f lush(stdout) ; 


} 

exit (0) ; 
#endif 


return  0; 

> 
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C.15  Simulation  of  the  Accumulator 


The  high-pass  filtered  image,  y,  provides  input  to  an  edge  detection  process  simulating  the  simple 
cells  of  area  VI.  Each  simple  cell  produces  a  number  of  outputs,  one  for  each  orientation  sensitivity. 
Thus  for  each  point  in  the  image,  a  vector  of  responses  is  produced,  and  each  position/orientation 
pair  simulates  a  single  axon  to  a  single  neuron  in  the  accumulator.  The  image  denoted  as  V  in  the 
flow  chart  is  therefore  a  vector  valued  function: 

v(r,  c) 

Even  though  we  now  have  a  three-dimensional  data  structure,  it  is  convenient  to  continue  to  call 
this  an  “image,”  and  the  elements  “pixels.” 

In  the  biological  accumulator,  all  SI  cells  are  simultaneously  active,  and  therefore  all  cells  in  the 
accumulator  may  receive  concurrent  stimulation.  In  computer  simulation  however,  we  are  forced 
to  scan  over  v.  The  simulation  process  is  accomplished  by  the  following  algorithm: 

for  (row  =  0;row<  nr  ;  row++) 
for  (  col  — 0;  col<  nc;col+- 1-) 

for(theta  =  0;theta  <  k ; theta=theta+dtheta ) 

{ 

t l=Gabor  (row  ,col  ,  theta) 
p  =  pointerto  (row  ,  col  ,  theta) 

d(rho  ,  theta)  =  Increment AccumulatorNeighborhood  (p  ,  tl ) 

{ 

The  pointers  are  computed  in  the  initialization  phase  of  the  simulation.  Each  pointer  simulates  a 
single  axon  from  the  SI  cell  at  (row,  col)  with  orientation  sensitivity  theta.  In  a  biological  neuron, 
the  accumulation  process  occurs  as  a  result  of  transfer  of  extracellular  sodium  into  the  cell  from 
numerous  synaptic  inputs.  Here  we  assume  the  presence  of  an  interneuron  which  accumulates  these 
inputs  and  produces  the  signal  d,  which  will,  in  turn  provide  stimulus  to  the  accumulator.  This  is 
simulated  by  simple  addition  in  the  function  IncrementAccumulatorNeighborhood,  which  adds  the 
signal  tl  to  the  interneuron  pointed  to  and,  in  a  Gaussian-weighted  way,  to  the  neighborhood  of 
that  interneuron. 

This  process  occurs  in  blocks  III-IV  of  the  flow  chart. 


/*=========================================================*/ 

/*  Accumulate  */ 

/*=========================================================*/ 


In  the  Accumulate  function,  the  function  is  passed  an  argumen  specifying  which  of  several  options 
to  use  for  accumulation 


int  Accumulate (IFSIMG  v,  float  **  Acc, struct  param  *p,int  mode) 
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{ 


double  dtheta; 
double  theta; 
int  i , j ; 

int  c;  //  loop  over  columns 
int  r;  //  loop  over  rows 

int  f;  //  loop  over  frames  (same  as  orientations 

int  nr;//  number  of  rows 

int  nc;  //  number  of  cols 

int  nf ; 

int  flatflag; 

int  negflag; 

int  irho;  //  integer  version  of  rho 

float  vtemp;  //place  to  store  a  temporatry 

double  sth,  cth;  //  place  to  remember  cos  and  sin  of  theta 

double  rho;  //the  rho  in  the  equation  of  a  st  line 

float  yl ,y2 ,y3 , a,b, cee ,xhat ,yhat ; 

int  xl,x2,x3; 

float  temp2; 

float  *mu; 

float  sum, average; 

double  Gx,Gy; 

double  mpl,mp2,spl; 

int  fmpl ,fmp2 , itheta; 

float  mul80[180]; 

void  FuzzAccumulator (struct  param  *) ; 

void  incacc (float  , float  **  ,  int  ,  int, int  .struct  param  *  , int); 
void  blurmu(float  *, float  *,  struct  param  *p) ; 
extern  void  Acc2ifs (float  **,char  * , int) ; 

mu= (float  *)malloc (p->numorientations  *  sizeof (float) ) ; 
nf  =  if sdimen(v,2) ;  //  v  is  a  3D  image 

nr  =  if sdimen(v, 1) ;  //  v  is  a  3D  image 

nc  =  if sdimen(v,0) ;  //  v  is  a  3D  image 

//  ifspot(p->v, "checkv. ifs") ; 
for(r=0;r  <  nr;r++) 

for(c=0;c<  nc;c++) 

{ 


The  following  code  estimates  uses  several  different  appraoches  to  incrementing  the  accumulator, 
given  the  18  (or  so)  measured  values  at  each  point  r,  c. 

First,  extract  all  18  measurements  and  find  the  average. 

sum  =  0.0; 

for (f =0 ; f <  p->numorientations ; f ++) 
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sum  +=  mu[f]  ; 

average  =  sum  /  p->numorientations ; 

if  the  Q  switch  is  defined,  rescale  the  values 

#define  Q 
#if def  Q 

f or (f =0 ; f <  p->numorientations ; f ++) 

{ 

mu[f]  =  mu[f]  *  mu[f]  /  average; 

} 

#endif 

Now  find  the  maximum  and  minimum 


mpl=  0 . 0 ; fmpl=0 ;mp2=0 ; fmp2=0 ; spl=  1000000.0; 
for (f =0 ; f <  p->numorientations ; f ++) 

{ 

mu[f]=  if sf gp3d(v,f , r , c) ;  //  get  Gabor  value  for  this  pixel 
if(mu[f]  >  mpl){mpl=mu[f] ;  fmpl  =  f;}  //remember  biggest 
if(mu[f]  <  spl)  spl=mu[f] ;  //  also  remember  smallest 

} 


#if def  DEBUG 

if (r  ==  c  &&  mpl  !=  0.0) 

{ 

printf("Mu  is  \n") ; 

for(f=0;f<p->numorientations;f++)  printf("°/of  ",mu[f]); 

} 

if  (r  ==  c&&  mpl  !=  0.0)  printf  ("mpl=°/0f  at  °/0d\n" , mpl , fmpl) ; 

#endif 

f or (f =0 ; f <  p->numorientations ; f ++) 

{ 

if (f  !=  fmpl  &&  mu[f]  >  mp2){mp2=mu[f] ;  fmp2  =  f;} 

} 

//  if(r==115  &&  c  ==  115)  printf  ("mpl  is  °/0f",mpl); 

now,  fmpl  is  the  angle  of  the  brightest  and  fmp2  is  the  second  brightest 

//  if (r  ==  115  &&  c  ==  115)printf ("jumping  to  mode  5\n"); 

switch  (mode) 

{ 

Version  0  is  the  classical  HT.  At  each  point,  the  entire  accumulator  is  incremented.  First,  the 
maximum  strength  of  the  response  is  tested. 
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Note  that  theta  ranges  from  -PAD  to  180+PAD  degrees.  This  is  because  the  accumulator  is 
deliberately  oversized  to  allow  smoothing. 

case  0: 

if (mpl  >  5.0) 

for(itheta  =  -PAD;  itheta<  180+PAD ; itheta++) 

{ 

incacc (mpl , Acc ,r , c , itheta,p, 1) ; 

} 

break; 


This  mode  uses  the  p-^numorientations  measurements  of  directional  derivitive  magnitude  and  then 
linearly  interpolates  them  over  a  total  of  180  1-degree  measurements.  Those  are  then  used  to  call 
incacc  the  interpolation  is  used  to  take  account  of  the  fact  that  we  have  only  p-inumorientations 
directional  filters,  not  180. 

case  1 : 

{ 

if (mpl  >  5.0) 

{ 

blurmu(&mu[0] ,&mul80  [0] ,p) ; 

f or (i=0 ; i<180 ; i++) 

incacc (mul80  [i] ,Acc,r,c,i,p,l) ; 
for(i=-PAD; i<0; i++) 

incacc (mul80 [180+i] ,Acc,r,c,i,p,l) ; 
for(i=180; i<180+PAD; i++) 

incacc (mul80 [180— i] ,Acc,r,c,i,p,l); 

} 

> 

break; 

mode  2  which  uses  the  B  matrix  to  make  a  mean-squared  estimate  gradient  direction,  using  all 
(p-inumorientation)  measurements 

case  2: 

{ 


Don’t  forget,  only  NUMORIENTATIONS  measuremnts  have  been  made  Now  multiply  by  the  Beta 
matrix  to  get  the  estimate  of  the  gradient 
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Gx=Gy=0 . 0 ; 

for(i=0;i  <  p->numorientations ; i++) 

{ 

Gx  +=  p->beta[l]  [i+1]  *  mu[i]  ; 

Gy  +=  p->beta[2]  [i+1]  *  mu[i]  ; 

} 

Now  that  we  have  the  gradient,  its  angle  is  from  the  arctangent.  The  arctan  will  return  numbers 
between  -pi  and  pi  and  so  we  have  to  convert  to  degrees  between  zero  and  . 


theta=atan2(Gy,Gx) ; 

We  are  operating  in  only  the  range  0  to  180  degrees,  so  if  we  add  n  to  an  angle,  we  must  change 
the  sign  of  rho.  That  is  accomplished  by  setting  negflag  before  calling  incacc. 

if (theta  <  0)  {theta  +=  M_PI;  negflag  =  -1.0;}else  negflag  =  1.0; 
theta  *=  RAD2DEG;  //  convert  to  floating  version  in  degrees 
itheta  =  theta+0.5;  //integerize  and  round 

increment  acc  using  value  from  Gabor  at  this  angle 

//  printf  ("Accumulate2 :  found  an  angle  of  °/0d\n" ,  itheta)  ; 

incacc (mu [itheta/p->degreesper sample] , Acc , r , c , itheta , p , negflag) ; 
}//end  of  if  statement 
break; 


Mode  3  uses  only  a  single  direction,  the  direction  of  the  maximum  responding  single  cell 


case  3: 

{ 

if(mpl>  5.0) 

{ 

for(i=0;i  <  p->degreespersample ; i++) 

incacc (mpl , Acc ,r , c ,fmpl  *  (p->degreespersample)  +i,p,1.0); 

> 

} 

break; 

Mode  4  takes  the  (p-^numorientations)  mu  values  and  increments  the  accumultor  at  a  point  corre- 
sponging  to  each  of  by  them  an  amount  proportional  to  its  strength. 

case  4: 

if (mpl>20 . 0) 

{ 

for(i=0;  i  <  p->numorientations ; i++) 
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for(j  =  0; j <p->degreespersample ; j++) 

{ 

itheta  =  i  *  p->degreespersample  +  j ; 
incacc (mu[i] , Acc ,r , c , itheta,p, 1.0); 

} 

} 

break; 

Mode  5  finds  the  brightest  direction  and  estimates  the  direction  as  the  weighted  average  of  the  max 
and  its  two  nearest  neighbors 


case  5: 

if (mpl  >  10.0  &&  (mpl-  spl)>0.05) 

{ 


We  already  know  the  biggest  element  of  mu.  That  is  at  fmp  with  magnitude  nip.  So  now  look  for 
the  neighboring  one 

//  if (r  ==  115  &&  c  ==  115) 

//  { 

//  printf  ("Accumulateafterif  :ml=°/„f ,  fmp  =  °/od\n"  ,mpl  ,fmpl) ;  ff  lush(stdout)  ; 

//  exit (0) ; 

//  > 

if(fmpl  ==  0)  xl  =  p->numorientations-l ; else  xl  =  fmpl-1; 
if(fmpl  ==  p->numorientations-l)  x3  =0;else  x3  =  fmpl+1; 
x2=fmpl ; 

yl=mu[xl] ; y2=mpl ; y3=mu [x3] ; 
if((yl==y2)  &&  (y2==y3))  flatflag  =  1; 
else  flatflag=0; 

Now,  we  fit  a  parabola  to  these  three  points.  Let  x2=0,  then  xl  =  -r  x3=r  the  three  quadratic 
equations  become 


y  1  =  r2a  —  rb  +  c 

(C.4) 

to 

II 

o> 

(C.5) 

y3  =  r2a  +  rb+  c 

if  there  are  r  degrees/sample  with  solution  shown  below 

(C.6) 

if  (flatflag  ==  0) 

{ 

cee=y2 ; 

b=(y3-yl)/(2.0*  p->degreespersample) ; 

a=(y3  -  (p->degreespersample) *b-cee) / (p->degreespersample  *  p->degreespersi 
xhat  =  -b/(2.0  *  a); 

yhat  =  a  *  xhat  *  xhat  +b  *  xhat  +cee; 
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case 


// 

// 


// 


// 

// 

// 


// 

// 


if (xhat  <  0.0  &&  fmpl  ==  0  ) 
xhat  =  xhat+180.0; 

else 

xhat  =  xhat  +  p->degreespersample*  fmpl; 
itheta  =  xhat  +0.5; 

> 

else 

{ 

itheta  =  p->degreespersample  *  x2; 
yhat  =  0.0; 

} 

incacc(yhat,Acc,r,c,itheta,p, 1) ; 

}//end  of  if  statement 


break; 

6  finds  the  brightest  two  values  and  interpolates  between  them 

case  6: 

if (mpl  >  5.0) 

{ 

float  pi ,p2 ,ExpectedVal ; 
pl=mpl/ (mpl+mp2) ; 
p2=mp2/ (mpl+mp2) ; 

if(c==196)  printf("pl  =  °/„f  at  theta  =  7„d  p2=°/0f  at  theta=°/0d\n" , 
pl,fmpl,p2,fmp2) ; 

if (fmpl  ==  0  &&  fmp2==p->numorientations-l) 
fmp2  =  p->numorientations ; 
if (fmpl  ==  p->numorientations-l  &&  fmp2==0) 
fmpl  =  p->numorientations ; 

ExpectedVal  =  fmpl  *  pi  +  fmp2  *  p2; 

itheta  =  ExpectedVal  *  10.0; 
itheta  =  ExpectedVal  *  p->degreespersample ; 

if (r  ==  c) 

printf  ("Calling  incacc  with  °/0f,0/0d  °/0d  °/0d\n" , 
mpl, r,c, itheta) ; 
incacc (mpl , Acc ,r , c , itheta, p, 1)  ; 

> 

break; 

default : printf ("unrecognized  accumulation  mode\n"); 
exit (-1) ; 

}//end  switch 

}//  end  loops  over  r  and  c 
Acc2if s (p->Accl , "Beforefuzz . if s" , 2*p->Accrows) ; 

FuzzAccumulator (p) ; //  t 
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free  (mu) ; 

}//  end  of  Accumulate 
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Figure  C.2:  Definitions  of  angles.  The  bold  line  is  the  line  of  interest.  It  sits  at  a  normal  distance  p 
from  the  origin.  The  angle  6  is  measured  by  the  Gabor  function  simulating  the  orientation-sensitive 
components  of  SI  cells. 


C.16  incacc 

This  function  simulates  activating  a  single  neuron  by  the  output  of  a  single  v  cell  The  choice  of 
theta  has  already  been  determined.  Theta  measured  in  degrees  The  familiar  formulation  is  initially 
used  for  the  equation  of  a  straight  line, (see  figure  C.2). 

p  =  x  cos  7  +  y  sin  7 

where  p  is  the  distance  from  the  origin  along  a  normal  to  the  line,  and  7  is  the  angle  that  the 
normal  (not  the  line  itself)  makes  with  the  x  axis. 


/* - - - - - — - -=-=== - - - ==== - */ 

/*  incacc  */ 

/*===============================================================================*/ 


void  incacc (float  v, float  **  Acc,  int  r,  int  c,  int  igamma, struct  param  *p,int  negflag) 

{ 

int  irho ; 

float  rho , value , accvalue ; 

extern  double  mysin(int) ,  mycos(int); 
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Convert  7  to  degrees  and  integerize 

rho  =  c  *  mycos (igamma)  +  r  *  mysin(igamma) ; 
rho  =  rho  *  negflag; 
irho  =  (rho  +0.5); 

irho  +=  p->Accrows;  //since  rho  could  be  negative,  shift  it. 
irho  +=  PAD;  //  set  the  0  0  point  correctly 

//  if (v  >1.0  kk  (c  ==  r)) 

//  { 

//  printf  ("incacc : %f  °/od  °/0d  °/„d"  ,v,r,c, igamma)  ; 

//  printf  ("irho  =  °/0d  itheta  =  7od\n" ,  irho ,  igamma) ; 

//  > 

/*  now  irho  is  in  array  coodinates  */ 

/*  check  valid  values  */ 

if  (irho  <  -PAD  )  {printf  ("neg  rho  %d  °/od\n" ,  irho ,  igamma) ;  f  flush  (stdout) ;  exit  (-1) ; } 
if  (irho  >  2*p->Accrows  +  2  *  PAD)  {printf  ("incacc :  rho  too  big  Accrows  is  %d  °/0d  °/od\n"  ,p->Ac< 
if  (igamma  <  -PAD  ) {printf  ("gamma  neg  °/0d  °/0d\n" ,  irho ,  igamma)  ; fflush(stdout)  ;  exit  (-1) ;  } 
if  (  igamma>  (180+PAD) ) {printf  ("gamma  of  °/0d  is  too  big  °/0d\n" ,  igamma,  irho)  ; f f  lush(stdout)  ;  e: 
//  if (c==196  kk  v  >  26.0) 

//  printf  ("float  coordinates  are  °/0d  °/0d  %d  °/0d\n"  ,r,c,irho,igamma+PAD) ; 

accvalue  =  Acc  [irho]  [igamma  +  PAD] ; 

//  if(c==r)  printf ("incacc :  updating  accumulator  cell  %d  %d  %f  \n" , irho , (igamma+PAD) , accvali 
Acc [irho] [igamma+PAD]  =  v+accvalue; 

//  if(c==r)  printf  ("incacc :  updated  accumulator  cell  %d  °/0d  is  %f  \n" ,  irho ,  (igamma+PAD) ,  Acc  [: 

//  printf ("Leaving  incacc .."); fflush(stdout) ; 

> 
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C.17  Low-pass  filtering  the  Accumulator 


The  function  FuzzAccumulator  will  blur  the  accumulator.  The  result  will  be  stored  in  one  of  the 
accumulator  buffers,  and  then  pointers  will  be  swapped  to  avoid  copying. 

//====================================================*/ 

//  FuzzAccumulator  */ 

//  apply  a  low-pass  fileter  to  the  Acc  */ 

/*  result  will  be  stored  in  p->Acc2,  then  the  pointers*/ 

/*  swapped  */ 

/*  the  accumulators  are  arrays  of  floats,  so  they  have*/ 

/*  to  be  turned  into  ifs  images.  Fortunately,  that  just*/ 

/*  requries  some  pointer  manipulation  */ 

/*====================================================*/ 

void  FuzzAccumulator (struct  param  *p) 

{ 

int  r, c, nr, nc, nonmax; 
int  i , j ; 
int  inducks ; 

IFSIMG  hdrl,hdr2; 

float  **tmpimg; 

float  sum; 

int  len[3] ; 

int  passes; 

float  *alptr , *a2ptr ; 

extern  void  Acc2ifs (float  **,char  *,int); 

extern  void  saveacc2d (float  **, struct  param  *,int); 

//  printf ("Enteringing  fuzz\n" ) ; f f lush(stdout) ; 

//  blur  the  accumulator  Accl  into  Acc2  using  a  kxk  kernel 
for(r=PAD;r  <  2*p->Accrows  +  PAD;r++) 
for(c=PAD;  c  <  180  +PAD;c++) 

{ 

sum  =0.0; inducks=0 ; 

f or (i=-p->kernelradius ; i<=p->kernelradius ; i++) 

for  (  j=-p->kernelradius; j<=p->kernelradius; j++) 

{ 

sum  +=  p->Accl [r+i]  [c+j]  *  p->kernel [inducks++]  ; 

//  printf  ("Fuzz :  sum=°/of  \n" ,  sum) ;  f  flush  (stdout) ; 
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} 

p->Acc2[r]  [c] =sum/inducks ; 

> 

accumulator  2  is  the  blurred  accumulator.  Make  it  accl 

tmpimg  =  p->Accl ;p->Accl  =  p->Acc2;  p->Accl=tmpimg; 

Now  the  blurred  accumulator  is  in  Accl. 

Acc2if s (p->Acc2 , "BeforeBlur . if s" ,2*p->Accrows) ; 

Acc2if s (p->Accl , "Af terBlur . if s" ,  2*p->Accrows) ; 

//  saveacc2d(p->Acc2 ,p , 2) ; //  2  means  show  padding  too 

#if def  NONMAX 

//  run  nonmaximum  suppression  on  the  accumulator 
for(passes  =  0;passes<2;passes++) 
for(r  =  PAD;r<2*(p->Accrows)+PAD;r++) 
for(c=PAD;  c  <  180  +  PAD; C++) 

{ 

{ 

nonmax=0 ; 

for(i=-2; i<=2; i++) 

{ 

for(j=-2; j<=2; j++) 

{ 

if (p->Acc2 [r]  [c]  <  p->Acc2  [r+i] [c+j]) 

{ 

nonmax  =  1 ; 
break; 

> 

if (nonmax)  break; 

> 

} 

if (nonmax)  p->Acc2 [r] [c]  =  1.0; 

} 

> 


//  find  how  many  peaks  there  are 

for(r  =  PAD;r<2*(p->Accrows)+PAD;r++) 
for(c=PAD;  c  <  180  +  PAD; C++) 

{ 

if (p->Acc2  [r]  [c]  >=  3000.0) 

printf("Acc  of  °/„f  at  °/„d  (°/0d)  °/od\n"  ,p->Acc2  [r]  [c]  ,r, r-p->Accrows-PAD ,  c) ; f flush 

> 

#endif 

//  swap  the  images  Accl  and  Acc2 
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tmpimg  =  p->Accl ;p->Accl  =  p->Acc2 ;p->Acc2=tmpimg; 
//  printf ("Leaving  fuzz\n") ; f f lush(stdout) ; 

> 
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C.18  Making  a  Gaussian  blurring  kernel 


/*  */ 
/*===================================================*/ 

/*  makekernel  */ 

/*===================================================*/ 

/*  function  to  make  a  convolution  kernel  */ 

/*  */ 


float  makekernel (float  *kernel,int  kernelradius) 

{ 

float  sigma, sigmasq,rsq, sum; 
int  index; 
int  i , j ; 

index  =0 ; sigma=kernelradius ; sigmasq  =  sigma*sigma; 
sum=0 ; 

f or (i=-kernelradius ; i  <=kernelradius ; i++) 

f or ( j  =-kernelradius ; j  <=kernelradius ; j  ++) 

{ 

rsq  =  i*i  +  j  *j ; 

kernel [index++] =exp(-  rsq/sigmasq)  /  (2.0  *  M_PI  *  sigma); 
sum  +=  kernel [index-1]  ; 

} 

for(index=0; index< (2*kernelradius+l)  *(  2*kernelradius+l) ;  index++) 
kernel [index] =kernel [index]  /  sum; 

//  printf ("makekernel :  sums  to  7of\n",sum); 
return  sum; 

> 
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C.19  Set  an  image  to  all  zeros 


Because  the  program  may  run  for  a  long  time,  it  is  necesary  to  sometimes  set  an  ifs  image  to  all 
zeros. 


/*===================================================*/ 

/*  zap  */ 

/*===================================================*/ 

/*  function  to  set  an  */ 

/*  ifsimage  to  all  zeros  */ 


void  zap(IFSIMG  z) 

{ 

char  *ptr; 
register  int  i; 
int  numbytes; 

numbytes=if sdimen(z , 0)  *  if sdimen(z , 1)  *  z->ifsbpd; 
ptr=z->if sptr ; 
f or (i=numbytes-l ; i>=0; i — ) 

*ptr++  =  0.0; 


} 
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C.20  Initializing  an  Accumulator 


Since  an  accumulator  is  a  pointer  to  a  vector  of  pointers  to  floats,  they  cannot  be  zapped  with  the 
ifs  zap  funtion  above. 

/*===================================================*/ 

/*  zapacc  */ 

/*===================================================*/ 

/*  function  to  set  an  accumulator  */ 

/*  to  all  ones .  */ 

void  zapacc (float  **z, struct  param  *p) 

{ 

float  *fptr; 
register  int  i; 
int  numpixels ; 
int  nr,nc,r,c; 

numpixels=  (180  +  2*PAD)  *  (2*p->Accrows  +  2*PAD) ; 
fptr=&z [0] [0] ; 
f or (i=numpixels-l ; i>=0 ; i — ) 

*fptr++  =  1.0; 


> 
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C.21  Disk  reread 


Reads  an  image  from  disk  if  the  image  has  already  been  read  once. 


/*===================================================*/ 

/*  ifsreread  */ 

/*===================================================*/ 


/*  reads  a  file  from  disk  into  an  exisitng  image  */ 
/*  note:  THIS  ONLY  WORKS  IF  THE  IMAGE  IS  UNCOMPRESSED*/ 

int  if sReRead(char  *f ilename,  IFSIMG  img) 

{ 

//  int  filesize; 

//  FILE  *fp,  *f openO  ; 

int  nr, nc, size; 

IFSIMG  in; 


nr  =  if sdimen(img, 1) ;  nc  =  if sdimen(img, 0) ; 

//  size  =  nr  *  nc  *  img->ifsbpd; 

//  printf  ("reread:  %d  rows  °/0d  cols  equals  %d  bytes\n"  ,nr,nc,size) ; 
//  fp  =  fopen(f ilename, "rb") ; 

//  if (f p  ==  NULL) 

//  { 

//  printf ("Cannot  read  file  named  %s\n" , filename) ; 

//  exit (-1) ; 

//  } 

//  f seek(fp , 512 , SEEK_CUR) ; 

//  f read( (img+512) ,  size,  1,  fp) ; 

//  //  memcpy( (img+512) ,fp, size) ; 

//  //  img->if sptr=(char  *) (img  +  512); 

//  f close(fp) ; 

/*  read  the  new  input  image  */ 

in  =  if spin(f ilename) ; 

memcpy (img->if sptr , in->if sptr , size) ; 

if sf ree (in, IFS_FR_ALL) ; 

//  if spot (img, "new_img. if s") ; 

return  0; 

> 
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C.22  Writing  a  file 

Debuging  function  not  currently  used 


void  writef ile (IFSIMG  y,char  *  prefix, int  index) 

{ 

char  name [32] ; 

sprintf  (name ,  "°/0s°/0d.  if s"  .prefix,  index)  ; 
printf  ("wriitng  to  the  file  named  °/0s\n"  .name)  ; 
ifspot(y.name) ; 

> 
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C.23 


Test  for  image  null 


/*  */ 

/*  ifnullimage  */ 

int  ifnullimage (IFSIMG  f, float  value) 

{ 

int  nr,nc,numpixels,i; 
float  *ptr; 

ptr  =  (float  *)  f->ifsptr; 
nr  =if sdimen(f , 1) ; 
nc  =if sdimen(f , 0) ; 
numpixels  =  nr  *  nc; 

//  printf  ("ifnullimage :  °/0d  rows,  °/0d  columns,  °/0d  floats\n" , nr  ,nc .numpixels)  ; 
for(i=0;i  <  numpixels ; i++) 

if(*ptr++  !=  value)  return  0; 

return  1 ; 

> 
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C.24  Clip  a  floating  point  image 

converts  a  float  image  to  a  float  image  with  values  between  o  and  255  Not  used  in  this  version 


/  * - :=---= - ----= - := - *7 

/*  clip  */ 

/•  . - .  -  -  . -  .  .  ■  ■  / 


void  clip (IFSIMG  a) 

{ 


float  *ptr; 
int  i , n ; 

n=if sdimen(a, 0)  *  if sdimen(a, 1) ; 
ptr  =  (float  *)a->ifsptr; 
f or (i=n-l ; i>=0 ; i — ) 

{ 

*ptr=255.0  *  1.0 

/ 

(1.0  +  exp(-  0.05  *  *ptr) ) ; 
ptr++; 

} 
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C.25  Blurring  the  vector  of  angle  measurements 


/• . - .  -  -  -  .  .  ■  ■  / 

/*  blurmu  */ 

/  •  - . - .  .  ---  -  ■/ 


void  blurmu (float  *mu, float  *mul80 , struct  param  *p) 

{ 

int  i , left , right ; 
float  mul,mur; 

//  for(i=0; i<170; i++) 

for(i=0;  i  < (p->numorientations-l) * (p->degreespersample) ; i++) 

{ 

left  =  (i/(p->degreespersample) ) * (p->degreespersample) ; 
right  =  (left+10) ; 

mul=mu [lef t/p->degreesper sample] ;mur=mu [right/p->degreespersample] ; 

//  if (mul  >  100  &&  mur  >  100) 

//  printf  ("Blumu:  l=%d  r=°/„d  i=0/0d  mu=7of\n" , left, right, i,mu[i] ) ; f f lush(stdout) ; 

mul80[i]=0.1  *  (mul  * (right  -  i)  +  mur  *(i  -  left)); 

} 

mul  =  mu  [p->numorientations-l] ;mur  =  mu[0] ; 

f or (i=(p->numorientations-l) * (p->degreesper sample) ; i< 180 ; i++) 

{ 

mul80[i]  =  mul  *  (180-i)  +  mul  *  (i- ( (p->numorientations-l) *p->degreespersample) ) 
//  if  (mul80  [i]  >  10000 . 0)printf  ("blurmu :mul80  [°/„d]  =  °/„f  \n" ,  i ,mul80  [i] ) ; 

} 


> 
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C.26  Display  the  Inverse  Accumulator 

/*=====================showHoughinv===================*/ 

/*  first,  call  Houghinv  to  get  the  inverse  HT  of  the  */ 

/*  HT  we  just  found  */ 

/*  then  display  it  and  hold  awaiting  a  keystroke  */ 

int  showHoughinv( struct  param  *p) 

{ 

extern  Houghinvsub (float  **,IFSIMG,  float  threshold, struct  param  *) 
extern  VideoscaleAcc (float  **, float  **, struct  param  *) ; 
int  nr,nc; 
int  displayindex; 

int  if nullimage(IFSIMG, float) ; 

int  WriteToIFSDisplayWindow(int , IFSIMG, int , float , int , int) ; 


//  printf ("Entering  showHoughinv: \n") ; f f lush(stdout) ; 
#undef  TESTHOUGHINV 
#if def  TESTHOUGHINV 

zapacc(p->Accl ,p) ; 

p->Accl [0+p->Accrows+PAD]  [135+PAD]  =  500; 

Houghinvsub (p->Acc 1 , p->reconstructed , 500 . 0 , p) ; 


#endif 

zap(p->reconstructed) ;  //  prepare  the  output  image 

Houghinvsub (p->Acc 1 , p->reconstructed , 500 . 0 , p) ; 
if (if null image (p->reconstructed, 2 . 0) ) 

{printf ("reconsturcted  is  empty") ; exit (0) 

//  p->display2  =  CreateIFSDisplayWindow(reconstructed, 1 . 0) ; 

WriteToIFSDisplayWindow(p->display2 ,p->reconstructed,0 , 1 . 0,300,0) ; 
ifspot(p->reconstructed, "reconstructed. ifs") ; 

> 
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C.27  MakeFake  Function  only  for  Testing 


/• . - .  -  -  -  . - .  .  ■  ■  / 

/*  makefake  */ 

/•-  - . - .  . ====  -  •  / 


/*  creates  a  synthetic  image  to  test  */ 
void  makefake (char  *argvl) 

{ 

IFSIMG  farble; 
int  len[3] ,row; 

len[0] =  2 ; len [1] =256 ; len [2] =256 ; 
farble=ifscreate("f loat" , len, IFS_CR_ALL , 0) ; 
for (row  =  50; row  <  100;row++) 

{ 

ifsfpp(farble, row, row, 200.0) ;  //  image  with  a  single  diagonal  line 
if sfpp(farble , row, 100-row, 200 . 0) ;  //  image  with  a  single  diagonal  line 

> 

if spot (farble , argvl) ; 

> 

C.28  Accumulator  Operations  Functions 

// 

//  Accoperations . c 

// 

// 

//  Created  by  Wesley  Snyder  on  6/13/12. 

//  Copyright  (c)  2012  _ MyCompanyName _ .  All  rights  reserved. 

// 

#def ine  TRUE  1 
#def ine  FALSE  0 

C.28.1  Video  Scale  Accumulator 


/  *  -  -  - .  - .  . - .  -  ■  •  / 

/*  VidscaleAcc  */ 

/*  */ 
/•  ---------  .  -  -  ---  ■  -  -  .  / 


/*  accepts  one  accumulator  in  standared,  padded*/ 
/*  format,  and  video  -scales  it  */ 

VidscaleAcc(f loat  **a,  float  **b, struct  param  *p) 

{ 

int  r , c ; 

float  max, min, scale ; 
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max  =  0.0;min=100000.0; 
for(r  =  PAD;  r  <  2*p->Accrows+PAD;r++) 
for(c  =  PAD;  c  <  180+PAD;c++) 

{ 

if  (a[r]  [c]  >  max)  max  =  a[r]  [c]  ; 
if (a[r]  [c]  <  min)  min  =  a[r]  [c] ; 

> 

scale  =  255.0  /  (  max-min) ; 
for(r  =  PAD;  r  <  2*p->Accrows+PAD;r++) 
for(c  =  PAD;  c  <  180+PAD;c++) 

b[r][c]  =  scale* (a [r] [c] -min) ; 


C.28.2  Mark  POints  of  Interest  in  Accumulator 


/*=============================================*/ 

/*  markaccl  */ 

/*  multiply  the  point  in  the  accumulator  to  -1  */ 

/*  so  it  can  be  identified  later  */ 

/*  the  marking  is  done  in  second  acc  to  avoid  */ 

/*  in-place  difficulties  */ 

void  markacc2 (struct  param  *p,int  r,int  c) 

{ 

p->Acc2  [r+p->Accrows  +  PAD]  [c+PAD]  = 
p->Accl [r+p->Accrows  +  PAD]  [c+PAD]  *  -1.0; 

> 

void  remarkaccl (struct  param  *p) 

{ 


> 


int  r , c ; 

for(r=0;r  <  2*p->Accrows  +  2*  PAD;r++) 
for(c=0;c  <  180  +  2  *  PAD;c++) 
if (p->Acc2  [r]  [c]  <0.0) 

p->Accl[r]  [c]  =  p->Acc2  [r]  [c] ; 


C.28.3  ShowAccNeighborhood 

Shows  the  3x3  neighborhood  of  a  point 

void  ShowAccNeighborhood(struct  param  *p, float  **Acc,int  r,  int  c) 

{ 

int  i , j ; 

printf  ("\n***°/0d  °/0d***\n" ,  r ,  c)  ; 
for(i  =  — 1 ; i  <=  1  ;i++) 

{ 


f or ( j=-l ; j<=l ; j++) 

printf("°/0f  ",Acc[r  +  p->Accrows  +  PAD  +i]  [c+PAD+j] )  ; 
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printf ("\n") ; f f lush(stdout) ; 

} 

> 

C.28.4  Determine  if  a  point  is  a  Local  Maximum 


/*=============================================*/ 

/*  localmax  */ 


/*  returns  TRUE  if  the  point  of  interest  is  a  */ 

/*  local  maximm  */ 

/*  note  that  r  and  c  are  values  in  rho  and  theta  */ 

/*  so  that  r  ranges  from  -Accrows  to  +accrows  */ 

/*  this  function  adds  the  PAD  */ 

int  localmax (struct  param  *p, float  **Acc,int  r,  int  c) 

{ 

int  nonmax; 
int  i ,  j  ; 

float  centervalue; 
nonmax=0 ; 

// 

centervalue  =  Acc [r+p->Accrows+PAD]  [c+PAD]  ; 

/* 

f or (i=r+(p->Accrows)+PAD-l ; i<=r+(p->Accrows)+PAD+l ; i++) 
f or ( j=c+PAD-l ; j<=c+PAD+l ; j++) if ( ! (i==r  &&  j  ==c)) 
if  (centervalue  <=  Acc[i][j]) 

*/ 

f or (i=-l ; i<=+l ; i++) 

f or ( j=-l ; j<=+l ; j++) if ( ! (i==0  &&  j  ==0)) 

if (centervalue  <=  Acc  [i+r+p->Accrows  +PAD]  [j+c+PAD] ) 
return  FALSE; 

//  printf  ("localmaxacc  of  °/0f  at  %d  °/0d  ",  centervalue,  r,c)  ;  f  f  lush(stdout)  ; 

//  printf  ("which  is  %f  at  %d  °/0d\n" ,  centervalue , 

//  r+p->Accrows+P AD , c+PAD) ; f f lush(stdout) ; 

return  TRUE; 


} 

C.28.5  Make  an  Accumulator 

float  **  MakeAccumulator (int  rows)  //  accumulators  have  180  columsn  (+pad) 

{ 

float  *Acc; 
float  **ptr; 
float  *p; 
int  i ; 

void  *  malloc(size_t) ; 
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Acc  =  (float  *)malloc ( (2*rows  +  2  *  PAD)  *  (180  +  2  *  PAD)  *  sizeof (float) ) 
ptr  =  (float  **)malloc( (2*rows  +  2  *  PAD)  *  sizeof (float  *)); 
p=Acc ; 

for(i  =  0;  i  <  2*(rows  +  PAD) ; i++) 

{ 

ptr  [i]  =  p; 
p  +=  (180  +  2  *  PAD); 

} 

//  ford  =0;i<  180 ;  i++)ptr  [rows]  [i]  =  128.0; 

//  now  the  value  returned  by  MakeAccumulator  can  be  accessed  using  []  [] 

//  printf ("Returning  from  MakeAccumulatpor") ; f f lush(stdout) ; 
return  ptr; 


C.28.6  test  Indexing 

This  function  is  a  debugging  function. 


/*=============testaccindexing=================*/ 

void  testaccindexing (struct  param  *p) 

{ 


int  rows , columns ; 
int  r , c ; 
float  v; 

for(r  =  0;  r  <  p->Accrows;r++) 
for(c  =  0;  c  <  180;  C++) 

v=p->Accl [r+PAD]  [c+PAD] ; 

> 

/*==========ReadAccumulator  =============*/ 

/*  coordinates  range  (in  degreess)  from  -pad  to  180  +  pad*/ 


float  ReadAccumulator (float  **A, float  rho,  int  itheta, int  rhos) 

{ 

int  irho ; 
void  exit (int) ; 

irho  =  rho  +0.5;  //  round  off  rho 

irho  +=  rhos  +  PAD; 


/*  angles  are  also  in  degrees,  ranging  from  minus  PAD  to  180+PAD*/ 


itheta  +=PAD; 

if(itheta  <  0  )  {printf  ("RD:  Invalid  Acc  address  70d\n" ,  itheta)  ;  exit  (-1) ; } 
return  A [irho] [itheta] ; 

> 

/*==========WriteAccumulator  =============*/ 

/*  coordinates  range  (in  degreess)  from  -pad  to  180  +  pad*/ 
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float  WriteAccumulator (float  **A, float  rho,  float  theta, int  rhos, float  value) 

{ 

int  itheta,  irho; 

irho  =  rho  +0.5;  //  round  off  rho 

irho  +=  rhos  +  PAD; 

/*  angles  are  also  in  degrees,  ranging  from  minus  PAD  to  180+PAD*/ 
itheta  =  theta  +  PAD; 


return  A [irho] [itheta]  =  value; 

> 

C.28.7  Cosine  of  an  Angle  Specified  in  Degrees 

/*========================myCOS===========================*/ 

double  mycos(int  itheta) 

/*  finds  the  cosine  and  sine  of  an  integer  described  in  degrees*/ 
/*  this  will  later  be  done  by  lookup  table*/ 

{ 

double  t; 

t=itheta  /  RAD2DEG; 
return  cos(t); 

> 


C.28.8  Since  of  an  Angle  Specified  in  Degrees 

/*========================myCOs===========================*/ 

double  mysin(int  itheta) 

/*  finds  the  cosine  and  sine  of  an  integer  described  in  degrees*/ 
/*  this  will  later  be  done  by  lookup  table*/ 

{ 

double  t; 

t=itheta  /  RAD2DEG; 
return  sin(t) ; 

> 

C.28.9  Convert  an  Accumulator  to  an  IFS  Image 


/  *==  ===== .  ===  ■  .  ==*/ 

/*  */ 
/*  Acc2ifs  */ 

/  •  ==  . - . - . ------- .  ===  -  ■  / 
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/*  writes  an  accumulator  out  to  an  ifs  image  */ 
/*  Acc2if s(floataccumulator, filename  */ 
/*  usually  used  for  debugging  */ 
/*  display  includes  the  pad  */ 


void  Acc2ifs (float  **acc,  char  *f ilename , int  rows) 

{ 

IFSIMG  outimg; 
int  len[3] ; 
int  r , c ; 
float  *fptr; 
len[0]  =2; 

len [1] =180  +  2  *  PAD; 
len[2]=rows+  2*PAD; 

out img=if screate ( "float" , len, IFS_CR_ALL , 0) ; 

fptr= (float  *)outimg->ifsptr; 
for(r  =  0;r<len[2] ;r++) 
for(c=0;c<len[l] ; C++) 

*f ptr++  =  acc  [r]  [c] ; 
if spot (outimg, filename) ; 

> 

/ *================saVeacc2d===============*/ 

/*  if  padflag  ==  1,  make  an  image  showing  */ 

/*  the  padding,  if  padflag  ==  2  dont  */ 

/*  otherwise  error  message  */ 

void  saveacc2d(f loat  **acc, struct  param  *p,int  padflag) 

{ 

IFSIMG  temp; 
int  len [3] ,r,c; 
float  value; 
switch (padflag) 

{ 

case  1 : 

len [0]  =2 ; len [1] =180 ; len [2] =2*p->Accrows ; 
temp  =  ifscreate("f loat" , len, IFS_CR_ALL , 0) ; 
for(r=0;r<  2*p->Accrows ; r++) 
f or (c=0 ; c<180 ; C++) 

{ 

value=acc [r+PAD]  [c+PAD] ; 
ifsfpp(temp,r,c, value) ; 

} 

break; 
case  2: 

len [0] =2 ; len [1] =180+2*PAD ; len [2] =2*p->Accrows+2*PAD ; 
temp  =  ifscreate("f loat" , len, IFS_CR_ALL , 0) ; 
for (r=0 ; r<  2*p->Accrows+2*PAD ; r++) 
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f or (c=0 ; c<180+2*PAD ; C++) 

{ 

value=acc [r]  [c]  ; 
ifsfpp(temp,r,c, value) ; 

} 

break; 
default : 

printf  ("Illegal  pad  value  of  °/0d  passed  to  saveacc2d" , 
padf lag) ; 
exit (-1) ; 

}//  end  switch 

ifspot(temp, "savedacc2d. ifs") ; 
ifsfree(temp,IFS_FR_ALL) ; 

> 


Utility  function  disc.  Draw  a  disc  on  an  ifs  image 


/•  .  - . - .  •/ 

/*  disc  */ 

/•  -----  -  -  -  . -  -  -  -  -  •/ 


void  disc(IFSIMG  img.int  r,int  c,int  rad,int  brightness) 

{ 

float  dtheta,  theta, x,y, radius ,f rad; 
int  ix,iy; 
frad  =  rad; 

for(radius  =  2.0;radius  <=  frad;radius  +=  1.0) 

{ 

dtheta  =  0.5  /  radius; 

for(theta  =  0.0;  theta  <  (2.0)  *  (M_PI);theta  +=  dtheta) 

{ 

x  =  c  +  radius  *  cos  (theta) ; 
ix  =  x+ . 5 ; 

y  =  r  +  radius  *  sin  (theta) ; 
iy  =  y+ . 5 ; 

if sipp(img, iy, ix, brightness) ; 

//  printf  ("r=°/od  c=°/0d  theta  =  °/0f  x  =  %d  y=  %d  br=°/„d\n"  ,r,c, theta, ix,iy, brightness) 

} 

} 


} 

/*===================================================*/ 

/*  copyandconvert  */ 
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/*===================================================*/ 

/*  converts  a  float  image  to  a  uchar  image  */ 

/*  copyandconvert(floatimage,ubyteimeage  */ 

/*  used  in  displaying  accumulator  */ 


void  copyandconvert (float  **  imgl , IFSIMG  img2) 

{ 

unsigned  char  *uptr; 
float  *fptr; 
int  nr,nc; 
int  row, col, iv; 
float  max, min; 
float  vf,vi, scale; 

nr  =  if sdimen(img2, 1) ;nc  =  if sdimen(img2,0) ; 

//  printf  ("CopyConvert :  nr=0/0d  nc  =  °/0d\n"  ,nr  ,nc) ; 

//  uptr= (unsigned  char  *) img2->if sptr ; 

Find  max  and  min  of  Accumulator 

max  =  -1000000 . 0 ;min=1000000 . 0 ; 
for(row=0;row  <  nr;row++) 

for(col  =  0;col  <  nc;col++) 

{ 

if (imgl  [row]  [col]  >  max)  max  =  imgl [row]  [col] ; 
if (imgl [row]  [col]  <  min)  min  =  imgl [row] [col] ; 

} 

scale=255.0/ (max-min) ; 
for(row=0;row  <  nr;row++) 

for(col  =  0;col  <  nc;col++) 

{ 

vf =imgl [row]  [col] ; 

if (isnan(vf) )printf ("copyandconvert  nan  °/0d  °/0d  \n", row, col) 

vi=(vf-min)  *scale; 
if (vf  >  max  /  2.0) 

{ 

//  printf  ("peaks  at  %d  °/„d  °/0f->°/„f\n"  ,row,col,vf  ,vi) ; 

disc(img2,row,col,5,vi) ; 

} 

/ /  *uptr++  =  vi ; 

} 

//  printf ("writing  test  image  to  foo.ifs\n"); 
if spot (img2 , "f oo . if s") ; 

now  that  all  the  points  have  been  rescaled,  mark  the  biggest  with  discs 

//  uptr= (unsigned  char  *) img2->if sptr ; 
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//  for(row=0;row  <  nr;row++) 

//  for(col  =  O;col  <  nc;col++) 

//  { 

//  iv=*uptr++; 

//  iv=if sigp(img2,row, col) ; 

//  if (iv  >  128) 

//  { 

//  printf  ("iv=°/0d  Drawing  disc  at  °/0d  70d\n" , iv,row, col) ; 

//  disc (img2 ,row, col , 5 , iv) ; 

//  > 

//  > 

> 

C.29  Reconstruct  and  Display  the  Original  Image 

// 

//  Houghinvsub . c 

// 

// 

//  Created  by  Wesley  Snyder  on  5/24/12. 

//  Copyright  (c)  2012  _ MyCompanyName _ .  All  rights  reserved. 

// 

Compatible  with  version  2.0  of  SLDS 

/*  ============================Houghinvsub=============================*/ 

/*  */ 

/*  program  to  read  an  output  from  a  hough  transform  and  drawlines  */ 

/*  in  an  image  */ 

/*  usage:  Houghinvsub (Houghimg,  outimg)  */ 

/*  note:  by  convention,  Houghimg  normally  has  180  columns,  */ 

/*  (one  column  per  degree)  */ 

/*  Parameters:  d  theta,  the  angle  in  degrees  per  column  value  */ 

/*  is  usually  1,  but  this  program  will  compute  it,  just  in  case  */ 

/*  the  input  does  not  have  180  columns  */ 

/*  it  will  still  assume  the  angles  range  from  0  to  180  degrees  */ 

/*  anr  will  be  p->Accrows.  */ 

/*  p->tempacc  is  a  FLOAT,  and  must  be  vidscaled  before  conferting  */ 

/*  outimg  for  display  */ 


int  Houghinvsub (float  **Houghin, IFSIMG  outimg, float  threshold, struct  param  *p) 

{ 
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int  nr,nc; 
int  icase; 
int  len[3] ; 

double  rho, theta, dtheta, st, ct; 

double  x,y,max; 

int  iy,ix,i; 

int  irho , itheta ; 

float  v,vin; 

int  anr , anc ; 

float  min,mostneg; 

int  mostnegtheta,mostnegrho , lastrho , lasttheta; 
anr=  p->Accrows ; 
anc  =  180; 

dtheta  =  M_PI  /  anc; 

nr  =  if sdimen(outimg, 1) ;  nc  =  if sdimen(outimg,0) ; 

First,  find  the  most  negative  point  in  the  accumulator  and  remember  it. 
lastrho  =  361; 

for(i=0;i  <  p->numpeaks ; i++) 

{ 

mostneg  =  0.0; 

for  (irho  =  -anr; irho  <  anr ; irho++)for(itheta=0;  itheta  <  anc ; itheta++) 

{ 

v=Houghin [irho+anr+PAD] [itheta+PAD]  ; 
if (v  <  mostneg  ) 

{ 

mostneg  =  Houghin [irho+anr+PAD] [itheta+PAD] ; 
mostnegrho  =  irho; 
mostnegtheta  =  itheta; 

} 

> 

printf  ("°/„d  °/0d\n" , mostnegrho , mostnegtheta)  ; 

p->peak[i]  =  mostneg ; p->ipeakrho [i] =mostnegrho ;p->ipeaktheta [i]  =  mostnegtheta; 

Houghin [mostnegrho+anr+P AD] [mostnegtheta+PAD]  *=  -1.0;  //  unmark  this  point 

> 

for(i=0;i  <  p->numpeaks ; i++) 

Houghin [p->ipeakrho  [i]  +  anr  +PAD]  [p->ipeaktheta [i] +PAD]  *=  -1.0; 

//  printf  ("Houghinvsub :  strongest  response  is  at  rho, theta  =  (°/od  °/0d)  =  (°/0d  ’/,d)\n",mostm 
//I  now  loop  over  the  list  of  peaks  drawing  lines  when  needed 
for  (i=0;i  <  p->numpeaks ; i++) 

{ 

//irho  is  coordinates  in  the  accumulator 
//  so  irho  is  always  positive,  must  convert  to  find  rho 
irho  =  p->ipeakrho [i] ; itheta  =  p->ipeaktheta[i] ; 
if (irho  ==  0)  rho  =  irho  +  drand48 O/1000.0;  else 
rho  =  irho; 

theta  =  itheta  *  dtheta; 
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st  =  sin(theta) ;  ct  =  cos(theta); 
vin=Houghin[irho+anr+PAD] [itheta+PAD]  ; 


if (vin  <0.0) 

{ 

//  printf  ("Houghinvsub:°/Od  °/od  %f  \n" ,  irho ,  itheta, vin) ; f  f  lush(stdout) ; 

if (fabs(st)<=  0.001) 

{ 

f or (iy=0 ; iy<nr ; iy++) 

{ 

if (vin  ==  mostneg) 

if sfpp(p->tempacc , iy , irho , -mostneg  *  2.0  ); 

else 

if sfpp(p->tempacc , iy , irho , -vin) ; 

> 

} 

else  if(fabs(ct)  <=  0.001) 

{ 

//  printf  ("horizontal  rho=°/0f  irho  =  °/0d,  itheta=°/0d  (theta  =  %f)  , 

for(ix=0; ix<nc; ix++) 

{ 

if (vin  ==  mostneg) 

if sfpp(p->tempacc, irho, ix, -mostneg  *  2.0  ); 

else 

if  sf pp (p->tempacc , irho , ix , -vin) ; 

} 

} 

else 

{ 

#def ine  ISINX(x)  (x>=0  &&  x<fnc)?l:0 
#def ine  ISINY(y)  (y>=0  &&  y<fnr)?l:0 

float  xO ,y0 ,xl ,yl , alpha, length, dalpha, fnc ,f nr ; 

int  icase; 

fnr  =  nr; fnc  =  nc; 

//  printf  ("\ndrawing  line  °/0d  °/„d  brightness  =  °/0f  \n" ,  irho ,  itheta, • 

//  compute  where  the  line  crosses  the  X  axis 
icase  =  0; 

if ( (ISINX( (rho/ct) ) )  &&  (ISINY( (rho/st) ) )  )  icase=icase  I  1; 

else 

if ( (ISINX( (rho/ct) ) )  &&  (ISINX( (rho-nr  *  st)/ct)))  icase=icase  |  2; 
else 

if ( (ISINX( (rho/ct) ) )  &&  (ISINX( (rho-nc  *  ct)/st)))  icase=icase  |  4; 
else 

if ( (ISINY( (rho/st) ) )  &&  (ISINX( (rho-nr  *  st)/ct)))  icase=icase  |  8 
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else 


if ( (ISINY( (rho/st) ) )  &&  (ISINY( (rho-nc  *  ct)/st)))  icase=icase 
else 

if ( (ISINX( ( (rho-  nr  *  st)/ct)))  &&  (ISINY( (rho-nc  *  ct)/st’ 

switch(icase) 

{ 

case  1:/*  printf("case  l:Line  crosses  x=0  and  y=0\n");*/ 
x0=0 ; yO=rho/st ; yl=0 ; xl=rho/ct ; 
break; 

case  2:  /*printf ("case  2: Line  crosses  y=0  and  y=nr\n");*/ 
y0=0;x0=rho/ct ;yl=nr ;xl=(rho-nr  *  st)/ct; 
break; 

case  4:  /*printf ("case  3:Line  crosses  y=0  and  x=nc\n");*/ 
y0=0;x0=rho/ct ;yl=nr ;xl=(rho-nr  *  st)/ct; 
break; 

case  8:  /*printf ("case  4: Line  crosses  x=0  and  y=nr\n");*/ 
x0=0;y0=rho/st ;yl=nr ;xl=(rho-nr  *  st)/ct; 
break; 

case  16:  /*printf ("case  5:Line  crosses  x=0  and  x=nc\n");*/ 
x0=0;y0=rho/st;xl=nc;yl=(rho  -  nc  *  ct)/st; 
break; 

case  32:  /*printf ("case  6:Line  crosses  x=nc  and  y=nr\n") ; */ 
yO=nr;xO  =  (rho  -  nr  *  st)/ct;xl=nc;yl= (rho-nc  *  ct)/st; 
break; 

default :printf ("invalid  icase  value  of  °/„d  ",icase); 
printf  ("rho=°/of ,  st=”/„f  ct  =  °/of\n"  ,rho,st,ct) ; 
exit (-1) ; 
break; 

} 

//  compute  where  the  line  crosses  the  y  axis 


//  printf  ("line  from  x=7of  y=°/0f  to  x=°/0f  y=70f\n,  brightness=7„f "  ,x< 

//I  the  followoing  code  calculates  all  the  points  along  a  line  which  is 
//I  specified  by  its  end  points.  If  $x_0  and  x_l,  x  \in  \Re~n$  are  the  vector  < 

//la  line  in  an  n-dimensional  space,  any  point  $x$  on  that  line  between  the  ei 

//I  $x=\alpha  x_0  +  (1-  \alpha  x_l$ 

length  =  sqrt ( (xO  -  xl)*(xO-xl)  +  (yO-yl)  *  (yO-yl)); 
dalpha  =  1.0/length; 

for (alpha=dalpha; alpha  <  (1 . 0-dalpha) ; alpha  +=  dalpha) 

{ 

x=alpha  *  xO  +  (1.0-alpha)  *  xl ; 

y=alpha  *  yO  +  (1.0-alpha)  *  yl ; 


127 


iy=y+0.5;ix  =  x+0.5; 
if (iy<0  )  iy=0 ; 
if(ix<0  )  ix=0; 
if(iy>nr-l  )  iy=nr-l; 
if(ix>nc-l  )  ix=nc-l; 
if(isnan(y))  printf ("Ynan") ; 

if (vin  ==  mostneg) 

if sfpp(p->tempacc , iy , ix, -mostneg  *  2.0  ); 

else 

if sfpp(p->tempacc , iy , ix, -vin) ; 


} 

} 

} 

if svidscale (p->tempacc , out img , &max , &min , 0) ; 

//  printf  ("Houghinvsub:vidscale  returned  %f  %f  \n" , max, min) 
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